VOLVER

Tema 6 Inferencia estadística

INFERENCIA ESTADÍSTICA

La inferencia estadística es el conjunto de técnicas estadísticas que nos permiten llegar a conclusiones sobre una población a partir de los datos de una muestra de dicha población. La estadística por tanto parte de los datos para extraer unas conclusiones, es un método deductivo. Para dar soporte teórico a la estadística existe el cálculo de probabilidades, que lo que proporciona es un método inductivo, es decir desde una regla es capaz de llegar a la misma conclusión. Tanto para la estadística como para la probabilidad, existen dos corrientes:

Frecuentista: Se basa en la repetición continuada de los hechos. Cuanto más se repita un experimento, la probabilidad de que ocurra un suceso será más regular y su frecuencia relativa se acercará mucho a la probabilidad del suceso.

USAREMOS ESTE MÉTODO CUANDO NO SEPAMOS ANTES LAS PROBABILIDADES DE CADA SUCESO, POR EJEMPLO QUE UNA MONEDA ESTÉ TRUCADA Y NO SEPAMOS CÓMO, SÓLO POREMOS SABER SU PARCIALIDAD REPITIENDO, POR ELLO SE LLAMA CONOCIMIENTO A POSTERIORI

La interpretación clásica, mayoritaria por lo menos hasta ahora, define la probabilidad en términos de experimentación. Si repites un experimento un número infinito de veces y compruebas que en 350 de cada 1.000 ocasiones se ha producido un determinado resultado, un frecuentista diría que la probabilidad de ese resultado es del 35%. Basándose en esta definición, un frecuentista afirma que es posible asociar a cada evento una probabilidad de obtener un valor VERDADERO del mismo.

La aproximación clásica se basa por lo tanto en estudiar la probabilidad real de las cosas, tratando de determinar hasta qué punto una medición realizada sobre un conjunto de experimentos se aproxima a la probabilidad real que subyace. Es por ello que un frecuentista definiría probabilidad como una expresión matemática que predice en qué medida es verosímil que ciertos eventos ocurran basándose en el patrón observado hasta este momento.

El enfoque frecuentista considera que los datos son aleatorios debido a la naturaleza del proceso que los genera y que cada vez que repitamos el experimento obtendremos un resultado diferente. De la misma forma cree que la hipótesis será cierta o falsa para el proceso estudiado pero que debido a la naturaleza aleatoria de los datos nuestro análisis puede señalarla como falsa cuando es cierta (falso negativo) o verdadera cuando es falsa (falso positivo). Por ello el enfoque frecuentista nos proporcionará, además de una respuesta a la pregunta de si la hipótesis es cierta o falsa, una probabilidad denominada p-valor. Este número es un indicador cuantas veces tendríamos que repetir la toma de datos y el análisis estadístico para obtener el resultado opuesto. Nótese que, al ser los datos aleatorios, siempre puede ser que acabe usando un conjunto que me apunte en la dirección contraria a la que realmente se comporta la población total.

Bayesiana: Hace un anális basado en la certeza que se tiene de que un hecho ocurrirá. Las observaciones se usan para actualizar o inferir la posibilidad de una hipótesis. Es una estrategia complementaria a la frecuentista, ya que con cada iteración añadiríamos un conocimiento obtenido en esa iteración para enriquecer el cálculo de la probabilidad.

USAREMOS ESTE MÉTODO CUANDO SEPAMOS A PRIORI LAS PROBABILIDADES DE CADA SUCESO, POR EJEMPLO AL LANZAR UNA MONEDA PERFECTA, POR ESO SE LLAMA CONOCIMIENTO A PRIORI,SABEMOS ANTES QUE LA PROBABILIDAD ES 50%

Por el contrario, la interpretación bayesiana se basa en un conocimiento limitado de las cosas. Afirma que sólo asocias una probabilidad a un evento porque hay incertidumbre sobre el mismo, es decir, porque no conoces todos los hechos. En realidad, un evento dado, o bien ocurrirá (probabilidad=100%) o bien no ocurrirá (probabilidad=0%). Cualquier otra cosa es una aproximación que hacemos del problema a partir de nuestro conocimiento incompleto del mismo. El enfoque bayesiano se basa por lo tanto en la idea de refinar predicciones a partir de nuevas evidencias. Un bayesiano definiría probabilidad como la expresión matemática que mide el nivel de conocimiento que tenemos para hacer una predicción. Por lo tanto, para un bayesiano, estríctamente hablando es incorrecto decir predigo que hay un 30% de probabilidades de que ocurra el evento P, sino que debería decir basándome en el conocimiento actual que tengo, tengo un 30% de certeza de que P ocurrirá.

Los frecuentistas piensan en los parámetros como unos valores fijos pero desconocidos. La estimación se basa en la elección de aquellos valores de los parámetros que maximizan la probabilidad de observar los datos.

Por su parte, los bayesianos interpretan los parámetros como variables aleatorias cuya distribución de probabilidad es estudiada en base al Teorema de Bayes. La idea es simple, un bayesiano ha de tener una distribución subjetiva de los parámetros antes de ver los datos (a priori) que modificará en función de los datos que haya observado para obtener una distribución a posteriori que resumirá todo el conocimiento del investigador sobre los parámetros de interés dados los datos y sus creencias a priori.

Los bayesianos le dan la vuelta al enfoque frecuentista. Para ellos los datos son fijos, no aleatorios porque ¿no se han obtenido ya? ¿Qué tienen de aleatorios los valores anotados en esas tablas que hemos ido generando? En cambio la hipótesis es para ellos aleatoria y puede ser o no verdad. Los bayesianos no determinan si la hipótesis es cierta o falsa, únicamente calculan la probabilidad (en sentido bayesiano) de que sea cierta o falsa.

Finalmente, sabemos que hay muchos resultados frecuentistas que pueden obtenerse desde una perspectiva bayesiana aunque la filosofía subyacente sea diferente. Mínimos cuadrados ordinarios (MCO) es el ejemplo más sencillo. Dicho estimador frecuentista coincide exactamente con la media de la distribución bayesiana bajo unas creencias concretas en el marco del modelo lineal normal. Existen en la literatura muchos otros ejemplos.

El objetivo de las técnicas de inferencia estadística es contrastar la hipótesis nula (Ho) y la hipótesis alternativa (H1), y en este sentido se pueden cometer dos tipos de error:

Error tipo I o alfa --> Rechazar la hipótesis nula siendo esta verdadera.
Error tipo II o beta --> Rechazar la hipótesis alternativa siendo esta cierta. La probabilidad complementaria del tipo II se llama poder de contraste (p-value) y es 1-beta.Este valor está totalmente impactado por el tamaño de la muestra, ya que el intervalo de confianza se va reduciendo cuanto mayor es el espacio muestral, por tanto más probable es rechazar Ho de manera efectiva.

La base de la inferencia estadistica es la probabilidad. Esta, asigna un número entre 0 y 1 que significa la posibilidad de que sea cierto. La probabilidad es el modelo de inferencia que se usa habitualmente para hacer estimaciones sobre fenómenos aparentemente aleatorios.

Teniendo en cuenta la conexión entre ambos enfoques en la práctica y la posibilidad de obtener los mismos resultados bajo ambos enfoques, los bayesianos argumentan que su interpretación del problema es siempre más intuitiva y natural. Recordemos que un bayesiano proporcionaría conclusiones del tipo: hay un 95% de probabilidad de que el parámetro esté entre 0,3 y 0,8, sin embargo, un frecuentista afirmaría: si generáramos 100 muestras aleatorias del mismo tamaño y repitiéramos la estimación 100 veces, en 95 de ellas el parámetro estimado estaría entre 0,3 y 0,8.

Los contrastes de hipótesis parecen por tanto más naturales en el marco bayesiano. Para un bayesiano convencido, un p-valor o un t-ratio es una herejía, ya que lo único que necesitamos para contrastar hipótesis es tener la distribución a posteriori de los parámetros. Salvo en contadas excepciones, un frecuentista basa su contraste de hipótesis en el análisis asintótico, es decir, en calcular p-valores de la distribución asintótica del estimador (no del parámetro verdadero, que es un valor fijo). Esta distribución, mayoritariamente normal gracias a numerosos teoremas centrales del límite, es la que el estimador tendría si el investigador tuviera muchos más datos de los que realmente tiene. ¿es ésta una buena forma de hacer inferencia en general? Por el contrario, el bayesiano se basa en la distribución de los parámetros dados sus datos, que puede en principio tener cualquier forma no gaussiana.

Al admitir información adicional, los métodos bayesianos reconocen que cada problema es distinto y promueven que el procedimiento de análisis se adapte al problema en cuestión, y no al revés; como consecuencia, tienden a ser más flexibles. Sin embargo, para poder hacer esto posible, la implementación de las técnicas bayesianas usualmente requiere de un esfuerzo computacional muy alto. La mayor parte de este esfuerzo se concentra en el cálculo de ciertas características de la distribución final del parámetro de interés. Por ejemplo, hay que integrar para pasar de una distribución conjunta a una colección de distribuciones marginales que sean útiles para hacer inferencias sobre los parámetros de interés. En la gran mayoría de los problemas las integrales requeridas no pueden resolverse analíticamente, por lo que es necesario contar con métodos numéricos eficientes que permitan calcular o aproximar integrales en varias dimensiones.

El enfoque bayesiano no sustituye el estudio de las frecuencias, lo enriquece con experiencias anteriores ya que puede ir refinando hacia atrás los resultados obtenidos en cada iteración. Un ejemplo extremo es que si en los últimos 4.500 millones de años en Cuenca ha salido el sol por el levante, suponemos que la probabilidad de que mañana vaya a volver a salir el sol es más alta con un estudio bayesiano de lo que sería con un análisis de estadística clásica, que no da ningún peso a la historia pasada. Pero cuando un médico pronostica le quedan unos seis meses de vida, querríamos creer que el galeno nos está facilitando esta información no sólo en base a un análisis de las frecuencias para esta dolencia, sino que además incorpora su experiencia de otros diagnósticos similares, y en especial nuestro propio historial clínico. Las conclusiones son afectadas por la premisa, pero si ésta es bien formulada, más que un defecto sería una ventaja

Afortunadamente, durante la segunda mitad del siglo XX se comenzaron a desarrollar técnicas numéricas flexibles y eficientes basadas en métodos de simulación estocástica. Esto, aunado al desarrollo de la tecnología que ha dado lugar a una mayor capacidad de procesamiento y de almacenamiento de los equipos de cómputo gracias a la virtualización, el Cloud Computing y el Big Data, está provocando un auge de los métodos bayesianos.

Probablidad es la cantidad de población que resume una aleatoriedad.Dos herramientas muy útiles para el cálculo de probabilidad son la distribución de probabilidad y la función de distribución

Esto es, ante la pregunta de si una hipótesis es correcta, los frecuentistas responden un sí o un no pero apostillan con qué frecuencia obtendrías el resultado opuesto mientras que los bayesianos no responden directamente y se limitan a indicar cómo de seguros están de uno u otro resultado.

Haremos usos de ambos métodos dependiendo de lo que necesitemos calcular. Usaremos por ejemplos métodos frecuentistas para el cálculo del error muestral y comparación entre muestras. Mientras que usaremos métodos bayesianos para el cálculo de especificidad y sensibilidad.

TEORÍA DE LA PROBABILIDAD

Dado un experimento aleatorio (por ejemplo lanzar un dado) su medida de probabilidad es la cantidad de dicha población que mide la aleatoriedad.

Está muy bien explicado en el siguiente artículo, del que he tomado algunas notas:

https://blogs.20minutos.es/mati-una-profesora-muy-particular/tag/probabilidad-condicionada/

Específicamente, la probabilidad toma un posible valor del experimento dónde:

  • Assigna un número entre 0 y 1.
  • La probabilidad de que algo ocurra seguro es 1.
  • La probabilidad de la unión de dos conjuntos de posibles resultados que no tienen nada en común (son mutuamente excluyentes, por ejemplo que salga un 2 o un 4 en un lanzamiento) es la suma de sus respectivas probabilidades.

Para asignar una probabilidad a un suceso, echamos mano de la regla de Laplace

laplace

laplace

En un caso más ámplio, como en el lanzamiento de dados

laplace

laplace

El matemático ruso Kolgomorov formalizó estas reglas de probabilidad.

  • La probabilidad de nada ocurra es 0
  • La probailidad de que sí ocurra es 1
  • La probabilidad de un hecho es 1 menos la probabilidad del hecho opuesto
kolgomorov

kolgomorov

  • La probabilidad de que sucedan al menos uno de dos (o más) cosas que no pueden ocurrir simultáneamente es la suma de sus probabilidades (UNIÓN DE CONJUNTOS)
kolgomorov

kolgomorov

  • Si un evento A implica la ocurrencia de un evento B, entonces la probabilidad de que ocurra A es menor que la probabilidad de que ocurra B (INCLUSIÓN DE CONJUNTOS)
qownnotes-media-GBLpWb

qownnotes-media-GBLpWb

  • Para cualquiera de dos eventos (o más) la probabilidad de que ocurran es la suma de sus probabilidades menos su intersección. (INTERSECCIÓN DE CONJUNTOS)
qownnotes-media-waZsrq

qownnotes-media-waZsrq

kolgomorov

kolgomorov

qownnotes-media-mTAmEg

qownnotes-media-mTAmEg

\[ \begin{eqnarray*} A_1 & = & \{\mbox{Person has sleep apnea}\} \\ A_2 & = & \{\mbox{Person has RLS}\} \end{eqnarray*} \] \[ \begin{eqnarray*} P(A_1 \cup A_2 ) & = & P(A_1) + P(A_2) - P(A_1 \cap A_2) \\ & = & 0.13 - \mbox{Probability of having both} \end{eqnarray*} \]

Una Variable aleatoria es el resultado numérico de un experimento. Las variables aleatorias tienen dos naturalezas, o discreta o contínua.

Un ejemplo de variable discreta es todo lo que sea contar un número de cosas(por ejemplo machos y hembras, o cualquier otra variable categórica), un ejemplo de continua es la medida de tu peso durante un mes.

Supongamos que deseamos conocer algún parámetro poblacional, por ejemplo la edad media de una población. En el contexto de la estadística clásica (frecuentista), podríamos extraer al azar múltiples y diferentes muestras representativas de forma sucesiva y determinar el estimador más apropiado en cada una de esas muestras. Para este caso, el estimador más apropiado para la media poblacional es la media muestral. Cada una de las medias de estas muestras tomarían un determinado valor diferente para cada muestra (o lo que es lo mismo, el estimador es tratado como una variable en sí mismo).En estadística, un estimador es un estadístico (esto es, una función de la muestra) usado para estimar un parámetro desconocido de la población.

Suponiendo la normalidad de la variable edad, los valores que irían tomando estas medias muestrales se distribuirían según una curva de Gauss (algunos de estos valores se darían con mayor frecuencia y otros con menor frecuencia). El área bajo dicha curva de distribución de frecuencias del estimador es la probabilidad total, es decir, la unidad.

DISTRIBUCIÓN DE PROBABILIDAD

Función de probabilidad (función de masa de probabilidad, PMF en ingles DISCRETOS)

Función que determina la probabilidad de cada evento posible, cuando los posibles valores son DISCRETOS: ejemplo:

  • valor p(1) -> 0.3
  • valor p(4) ->0.6
  • valor p(7)->0.1
PMF valores discretos

PMF valores discretos

Las distribuciones de probabilidad en base a su PMF más comunes son Bernoulli (valor p o p-1 muy útil para la prevalencia de una enfermedad), binomial (bernuilli repetido n veces) o poisson(contar el número de veces que se da un evento en el tiempo)

PMF valores discretos

PMF valores discretos

Función de densidad (PDF en ingles CONTINUO)

Función que determina la probabilidad de cada evento posible, cuando los posibles valores son CONTINUOS, por ejemplo probabilidad asociada a tiempo de espera en urgencias.

Ejemplo la proporción de llamadas que pueden ser atendidas en un día para dar citas en consultas externas sigue la siguiente curva:

\[ f(x) = \left\{\begin{array}{ll} 2 x & \mbox{ for } 1 > x > 0 \\ 0 & \mbox{ otherwise} \end{array} \right. \]

x <- c(-0.5, 0, 1, 1, 1.5); y <- c( 0, 0, 2, 0, 0)
plot(x, y, lwd = 3, frame = FALSE, type = "l")

Cuál será la probabilidad de atender el 75% al menos de las llamadas? Esta gráfica representa su PDF

Esta operación es muy sencilla ya que se puede calcular con trigonometría o através de la función beta, que es justo para distribuciones de probabiidad trianangular

1.5 * .75 / 2
## [1] 0.5625
pbeta(.75, 2, 1)
## [1] 0.5625

Por tanto la probabilidad de atender al menos el 75% de las llamadas es del 56%

Ya que calculamos el área bajo la curva, si elegimos un valor concreto su pdf nos dara 0,ya que el area de una recta es 0.

Area bajo curva de una recta es 0

Area bajo curva de una recta es 0

Las distribuciones de probabilidad en base a su PDF más comunes son Normal (la normal estandar es la que tiene e(x) = 0 y d(x) = 1)o de Gaus o exponencial (el tiempo entre dos eventos)

Los histogramas, pueden representar esta función de densidad en base a probabilidades o bien en base a número de ocurrencias de cada evento. EN nuestro caso, esto es fundamental ya que las probabilidades que vamos a manejar en el mundo de la la sanidad están todas basadas en la ocurrencia de un suceso, no en un cálculo a priori de sus probabilidades.

A la hora de representarlo, el histograma será mucho menos fiable para conocer su PDF ya que depende en gran medida del número de separaciones que hagamos, es más fiable utilizar diagrama de densidad kernel

http://www.statmethods.net/graphs/density.html

FUNCION DE DISTRIBUCION(CDF en ingles)

Es la probabilidad igual o menor asociada a un valor (probabilidad acumulada). Nos va a resultar muy útil para hacer el cálculo de la media en variables contínuas, ya que para su cálculo necesitamos este valor asociado a la variable.

\[ F(x) = P(X \leq x) \] Por otro lado la función de supervivencia será la que represente la probabilidad mayor a un valor

\[ S(x) = P(X > x) \]

Como se puede ver por tanto \(S(x) = 1 - F(x)\)

En variables contnuas, la PDF es la derivada de la CDF, como se puede ver en el ejemplo:

\[ F(x) = P(X \leq x) = \frac{1}{2} Base \times Height = \frac{1}{2} (x) \times (2 x) = x^2 \] La función de supervivencia por tanto será: \[ S(x) = 1 - x^2 \]

Esto puede ser calculado facilmente con la función pbeta, por ejmplo la probabilidad de atender el 40,50 y 60% de las llamadas

pbeta(c(0.4, 0.5, 0.6), 2, 1)
## [1] 0.16 0.25 0.36

Las probabilidades de la CDF se puede expresar también en Cuantiles y Percentiles

CUANTILES cuando tienes cuantil 0.95 significa que el 95% de la población está por debajo de tí, es decir se centra en el valor de la probabilidad menor o igual en lugar del resultado.

El cuantil alfa es el valor de la X que te da un valor concreto de cuantil. Ejemplo Cuantil alfa de 0.98 es X=6 PERCENTIL es el cuantil expresado en porcentaje

La media es el percentil 50%

En este ejemplo, queremos saber el cuantil 0.5 es decir la probabilidad de atender el 50% de las llamadas \(0.5 = F(x) = x^2\)

sqrt(0.5)
## [1] 0.7071068

Podemos hacer uso téambin de nuevo de pbeta:

qbeta(0.5, 2, 1)
## [1] 0.7071068

Hasta aqui hemos hablado en todo momento de medias sin tener en cuenta los datos, esto es así porque estamos calculando la media de la POBLACIÓN, un modelo de probabilidad conecta los datos con la población a través de asunciones. Por tanto la media de la que hemos estado hablando es el estimando, mientras que la media de la máuestra será el estimador.

PROBABILIDAD CONDICIONAL

La probabilidad de que un suceso pase sabiendo que ha pasado otro. En este caso, el hecho de conocer de antemano un hecho relaccionado (DEPENDIENTE) nos hace modificar la probabilidad de que ocurra el hecho en el que estamos interesados.

Esta es la base de la aproximación Bayesiana, ya que aportamos conocimiento a nuestra predicción y modificamos su probabilidad respecto a ese conocimiento.

  • Pongamos un hecho B = de manera que \(P(B) > 0\)
  • La probabilidad condicionada de otro evento A, sabiendo previamente que ha ocurrido B es definido como la probabilidad de que los DOS ocurran dividido por la probabilidad de que B occurra \[P(A\:|\:B) = \frac{P(A \:\cap\: B)}{P(B)}\]
  • En el caso de A y B sean independientes, entonces \[P(A\:|\:B) = \frac{P(A)P(B)}{P(B)} = P(A)\]

Un ejemplo de probabilidad condicionada, es qué probabilidades tengo de sacar un 1 si sé que ha salido impar.

qownnotes-media-jvqfcq

qownnotes-media-jvqfcq

El típico ejemplo de sucesos dependientes es el de sacar bolas de una bolsa. Si tenemos 5 bolas rojas y 3 bolas verdes en una bolsa, ¿cuál es la probabilidad de meter la mano y sin mirar sacar una bola roja?

sucesosDependientes

sucesosDependientes

La dejamos fuera ¿Y la probabilidad de que la segunda sea roja?Esto implica que el resultado del segundo experimento depende del primero. Lo que queremos conocer es la probabilidad de que en la segunda extracción (ya hay una bola fuera) saquemos una bola roja. A esa probabilidad que queremos calcular la llamaremos P(2ªR). Para calcular P(2ªR) vamos a necesitar calcular otras probabilidades antes. Por ejemplo, la probabilidad de que la 2ª bola sea roja condicionada al hecho de que la primera que salió fuese roja, lo llamamos P(2ªR | 1ªR)

sucesosDependientes

sucesosDependientes

Una vez tenemos todas las partes implicadas, será sumar las probabilidades de que la segunda haya sido roja condicionada a que la primera sea roja y que la segunda sea roja condicionando a que la primera no sea roja

sucesosDependientes Este caso es muy sencillo y podemos haer los cálculos símplemente contando

sucesosDependientes

sucesosDependientes

Ya sólo queda hacer la suma:

sucesosDependientes

sucesosDependientes

Teorema de Bayes Es el uso más famoso de la probabilidad condicional. Calcula la probabilidad de A condicionado a B desde la probabilidad de que pase B condicionado a A. Este teorema nos da la posibilidad de conocer las probabilidades condicionadas cruzadas de ambos eventos entre sí.

\[P(B\:|\:A) = \frac{P(A\:|\:B)P(B)}{P(A\:|\:B)P(B)+P(A\:|\:B^c)P(B^c)}\] donde \(B^c\) = corresponde a la probabilidad de que no suceda el evento \(B\), \(P(B^c) = 1 - P(B)\)

Es muy útil para calcular la efectividad de test diagnósticos ya que nos relacciona la efectividad del resultado de un test (Hecho A) con que realmente se tenga la enfermedad (Hecho B). Ejemplo: + y - son los resultados de diagnóstico de un test médico, y D que el sujeto tenga la enfermedad.

Probabilidad de que el test sea positivo = \(P(+)\) Probabilidad de que el test sea negativo = \(P(-)\) Probabilidad de que se tenga la enfermedad = \(P(D)\) Probabilidad de que no se tenga la enfermedad = \(P(D^c)\)

SENSIBILIDAD es la probabilidad de que el test sea positivo sabiendo que individuo tiene la enfermedad (Nos indica cómo ha funcionado). Es desable que sea alto, ya que indicará que nuestro test es fiable cuando da positivo.

Sensibilidad = \(P(+\:|\:D)\)

ESPECIFIDAD es la probabilidad de que el test sea negativo sabiendo que el individuo no tiene la enfermedad (Nos indica cómo ha funcionado). También es desable que sea alto, ya que indicará que nuestro test es fiable cuando da negativo.

Especificidad = \(P(-\:|\:D^c)\)

Prevalencia aparente Es la probabilidad de tener sepsis observada = \(P(D)\)

Prevalencia real es practicamente imposible de determinar pero se puede aproximar desde la prevalencia aparente con el estimador Rogan-Gladen

Valor predictivo positivo Probabilidad de que tenga la enfermedad sabiendo que el test es positivo (Nos indica cómo va a funcionar, tiene en cuenta prevalencia) = \(P(D\:|\:+)\)

Valor predictivo negativo Probabilidad de no tener la enfermedad sabiendo que el test es negativo (Nos indica cómo va a funcionar, tiene en cuenta prevalencia) = \(P(D^c\:|\:-)\)

Razón de verosimilitud diagnóstica -> Existen con el fin de dimensionar el beneficio clínico de un test de manera independiente a la prevalencia, ya que esta puede variar.

Razón de verosimilitud diagnóstica de test positivo El ratio entre la posibilidad de que en un paciente séptico haya saltado la alerta y que en un paciente no séptico haya saltado la alerta \[DLR_+ = \frac{sensitivity}{1-specificity} = \frac{P(+\:|\:D)}{P(+\:|\:D^c)}\] \(DLR_+ = .997 / (1 - .985) \approx 66\) DLR+ Ratio de 66 significa que la probabilidad de tener la enfermedad después del test es 66 veces mayor que antes del test

Razón de verosimilitud diagnóstica de test negativo Es el ratio entre la posibilidad de que en un paciente séptico no haya saltado la alerta y la posibilidad de que en un paciente no séptico no haya saltado la alerta. \[DLR_- = \frac{1 - sensitivity}{specificity} = \frac{P(-\:|\:D)}{P(-\:|\:D^c)}\] \(DLR_- = (1 - .997) / .985 \approx .003\) DLR- es 0.005 significa que la probabilidad de no tener la enfermedad antes del test es 0.003 veces mayor que después del test

Odds ratio Asociación entre exposición y resultado. Ejemplo de razón de oportunidad de indice de cancer de pulmón en fumador es 2. Esto significa que fumadores tienen 2 veces más probab. de tener cancer. = \(p/(1-p)\)

Diferencia entre probability y odds:

Probabilidad 0.80 Odds 0.80/1-0.80 = 4 –> Odds de 4 a 1

* $\frac{P(D)}{P(D^c)}$ = **pre-test odds**, o razón de oportunidad de enfermedad en ausencia del test
* $\frac{P(D\:|\:+)}{P(+\:|\:D^c)}$ = **post-test odds**, o razón de oportunidad de enfermedad dado un test positivo
* $DLR_+$ = factor por el cual la razón de oportunidad en la presencia de un test positivo puede ser multiplicado para obtener el post-test odds
* $DLR_-$ = relacciona la disminución en la razón de oportunidad de enfermedad después de un test negativo
* Siguiendo el ejemplo, para una sensibilidad de 0.997 y especificidad de 0.985, los razones de verosimilitud diagnosticos son $$DLR_+ = .997/(1-.985) = 66 ~~~~~~ DLR_- =(1-.997)/.985 = 0.003$$
* Esto indica que en el resultado de un test positivos, la razón de verosimilitud de la enfermedad es 66 veces el pretest odds
Test Séptico No séptico TOTAL
Pos 8503 2758 11261
Neg 37 355024 355061
TOTAL 8540 357782 366322

Prevalencia aparente= (8503+37)/(8503+37+2758+355024)=0.023

Prevalencia real = (Prev.apa+Espec-1)/(Sensib.+Espec-1) = (0.023+0.992-1)/(0.995+0.992-1)= 0.01519757

Proporción falso positivo = 2758/11261 = 0.24

Proporción falson negativo = 37/355061 = 0.0001

Sensibilidad = 8503/(8503+37) = 0.995

Especificidad = 355024/(2758+355024) = 0.992

VPP = 8503/(8503+2758) = 0.75

VPN = 355024/(2578+355024) = 0.992

DLR+ = (sensibilidad/1-especificidad) = 124

DLR- = (1-sensibilidad/especificidad) = 0,005

Odds ratio (razón de oportunidades) = (8503/2758)/(37/355024) = 29582.43

Ratio o razón –> Indice obtenido al dividir dos cantidades. Ninguno de los elementos del numerador estón incluidos en el denominador

Ejemplo Ratio alerta post = 8503/2758 = 3.08 por cada 3 alertas salta 1 que no es

Proporción -> El numerador está incluido en el denominador

    Proporcion alerta post =  8503/(8503+2758) = 0.75

Tasa -> Magnitud de cambio de un parámetro por unidad de cambio de otro que incluye una medida adiconal de tiempo en el denominador

Tasa de pacientes septicos en 2018: 8540/366322 = 0.023 es una tasa
de 23 pacientes septicos por cada 1000 urgencias.
qownnotes-media-tBAWRb

qownnotes-media-tBAWRb

POSITIVE PREDICTIVE VALUE es la probabilidad de tener la enfermedad sabiedo que el test es positivo.

NEGATIVE PREDICTIVE VALUE es la probabilidad de no tener la enfermedad sabiendo que el test es negativo.

En la ausencia de un test, se llama PREVALENCIA a la probabilidad de tener esa enfermedad.

qownnotes-media-dOUpkP

qownnotes-media-dOUpkP

qownnotes-media-tLsDNl

qownnotes-media-tLsDNl

qownnotes-media-pPYGcB

qownnotes-media-pPYGcB

qownnotes-media-RUIsxz

qownnotes-media-RUIsxz

DLR+ es cada cuantos test positivos acertados se da un test negativo err?neo en este caso 998/15 = 66, cada 66 test positivos acertados se da 1 negativo err?neo

DRL- es cada cuantos test positivos erróneos se da un test negativo acertado.

qownnotes-media-PaNiRC

qownnotes-media-PaNiRC

qownnotes-media-xTKrgI

qownnotes-media-xTKrgI

qownnotes-media-CtTfPs

qownnotes-media-CtTfPs

qownnotes-media-DynpEw

qownnotes-media-DynpEw

qownnotes-media-htWhsk

qownnotes-media-htWhsk

qownnotes-media-EtEyvi

qownnotes-media-EtEyvi

ESTUDIO DETALLADO DE LA SENSIBILIDAD Y ESPECIFICIDAD

Una vez que la prueba se ha aplicado a una población determinada, varios índices se utilizan para evaluar la prueba: - Verdadero positivos (VP): Número de casos que la prueba declara positivos y que son verdaderamente positivos. - Falsos positivos (FP): Número de casos que la prueba declara positivos y que en realidad son negativos. - Verdadero negativo (VN): Número de casos que la prueba declara negativos y que son realmente negativos. - Falsos negativos (FN): Número de casos que la prueba declara negativos y que en realidad son positivos. - N = TP + FP + FN + TN

qownnotes-media-HqjvcR

qownnotes-media-HqjvcR

  • Sensibilidad (equivalente a la tasa de positivos verdaderos): Proporción de casos positivos que están bien detectadas por la prueba. La definición matemática es: Sensibilidad = VP / (VP + FN).

  • Especificidad (también llamada Tasa de verdaderos negativos): proporción de casos negativos que son bien detectadas por la prueba. La definición matemática es: Especificidad = VN / (VN + FP).

  • Tasa de falsos positivos (FPR): Proporción de casos negativos que la prueba detecta como positivos.

  • Falso Negativo Precio (FNR): Proporción de casos positivos que la prueba detecta como negativo. - Prevalencia: la frecuencia relativa de los acontecimiento de inter?s en la muestra total (VP + FN) / N.

  • Valor Predictivo Positivo (VPP): Proporción de casos verdaderamente positivos entre los casos positivos detectados por la prueba. Tenemos PPV = TP / (TP + FP), o PPV = Sensibilidad x Prevalencia / [(Sensibilidad x Prevalencia + (1-Especificidad) (1-Prevalencia)]. Es un valor fundamental que depende de la prevalencia, un índice que es independiente de la calidad de la prueba.

  • Valor predictivo negativo (VPN): Proporción de casos verdaderamente negativos entre los casos negativos detectados por la prueba. Tenemos VPN = VN / (VN + FN), o VPP = Especificidad x (1 - Prevalencia) / [(Especificidad (1-Prevalencia) + (1-sensibilidad) x Prevalencia] . Este ?ndice depende tambi?n de la prevalencia que se independiente de la calidad de la prueba.

  • Razón de verosimilitud positiva (LR +): Esta relación indica a qué punto una persona tiene más posibilidades de ser positivo en la realidad cuando la prueba es que está diciendo es positivo. Tenemos LR + = sensibilidad / (1-especificidad). La RP + es un valor positivo o nulo.

  • Razón de verosimilitud Negativa (LR-): Esta relación indica a qué punto una persona tiene más posibilidades de ser negativo, en realidad, cuando la prueba es que está diciendo es positivo. Tenemos LR-= (1-sensibilidad) / (especificidad). El LR-es un valor positivo o nulo.

  • Odds ratio: El odds ratio indica la cantidad de un individuo es más probable que sea positivo si el resultado es positivo, en comparación con los casos en que la prueba es negativa. Por ejemplo, un odds ratio de 2 significa que la probabilidad de que el caso positivo se produce es dos veces superior si la prueba es positiva que si es negativo. El odds ratio es un valor positivo o nulo. Tenemos Odds ratio = VPxVN / (FPxFN).

  • El riesgo relativo: El riesgo relativo es un ratio que mide mejor la prueba se comporta cuando se trata de un informe positivo que cuando es negativo. Por ejemplo, un riesgo relativo de 2 significa que la prueba es dos veces más potente cuando es positiva que cuando es negativo. Un valor cercano a 1 corresponde a un caso de independencia entre las filas y columnas, y una prueba de que funciona tan bien cuando es positiva que cuando es negativo. El riesgo relativo es un valor nulo o un valor positivo dado por: riesgo relativo = VP / (VP + FP) / (FN / (FN + VN)).

VARIABLES INDEPENDIENTES

En este apartado cambiamos totalmente de premisa, ya que en los anteriores apartados las variables estaban relaccionadas entre sí, en este caso que ocurra la primera no infiere en que ocurra la segunda.

Dos variables son independientes si la ocurrencia de una no tiene nada que ver en la ocurrencia de otra.

Por ejemplo tirar dos veces una moneda y que salga las dos cara, son sucesos independientes 0.5x0.5 = 0.25

Si los sucesos no son independientes la intersección de las probabilidades no se pueden calcular multiplicando.

  • dos eventos \(A\) y \(B\) son independientes si lo siguiente se cumple:
    • \(P(A\:\cap\:B) = P(A)P(B)\)
    • De manera equivalente \(P(A\:|\:B) = P(A)\)
  • dos variables \(X\) e \(Y\) son independientes, si para cualquier subconjunto, A y B, se cumple
    • \(P([X \in A]\cap[Y \in B]) = P(X \in A)P(Y \in B)\)
  • independence = statistically unrelated from one another
  • si \(A\) es independiente de \(B\), entonces se cumple
    • \(A^c\) es independiente de \(B\)
    • \(A\) es independiente de \(B^c\)
    • \(A^c\) es independiente de \(B^c\)

Son sucesos independientes, por ejemplo, el lanzamiento consecutivo de monedas. Si empezamos en el punto azul y caminamos sobre las aristas azules, tenemos todos los resultados

SucesosIndependientes

SucesosIndependientes

Si queremos saber qué probabilidad hay de sacar 4 caras seguidas, no tenemos más que seguir el camino

SucesosIndependientes

SucesosIndependientes

Se trata de calcular la probabilidad de una intersección de sucesos independientes, porque el resultado de cada lanzamiento es independiente del anterior. Un suceso es que la primera salga cara, otro que la segunda salga cara, el tercero que salga cara en la tercera y el cuarto, que la cuarta tirada dé cara. Cada uno de estos sucesos, tiene probabilidad 1/2, como hemos dicho, Eso significa que en cada tirada de las 4, la probabilidad de que salga cara es 1/2 o sea, de un 50%.

SucesosIndependientes

SucesosIndependientes

Ahora bien, la probabilidad de que las cuatro veces nos salga cara es la probabilidad de la intersección, y eso se calcula sabiendo que, como son independientes, la probabilidad de la intersección es igual al producto de las probabilidades.

SucesosIndependientes

SucesosIndependientes

Es muy fácil caer aquí en la falacia del apostador o de Montecarlo. Este nombre viene porque en su casino, sucedió un hecho insólito hasta la fecha, una noche en la ruleta salieron 20 números negros consecutivos, por lo que las apuestas por roja en la siguiente tirada fueron en masa, con la creencia de que si habían salido 20 resultados seguidos negros, el siguiente sería muy difícil que fuera negro. Esa noche llegaron a salir hasta 26 veces seguidas números negros. Es la falsa creencia de que lo sucedido anteriormente en un experimento aleatorio, al azar, afectará al resultado de los experimentos siguientes.

Se están confundiendo la probabilidad de obtener 5 caras seguidas con la probabilidad de obtener cara en el 5º lanzamiento.

SucesosIndependientes

SucesosIndependientes

Todos los estudios que vamos a realizar a continuación, no sólo están basados en sucesos independientes, si no que además están identicamente distribuidos, IID variables

  • se dice que dos variables aleatorias son IID si son independientes e identicamente distribuidas
    • independientes = estadisticamente no relaccionadas entre sí
    • identicamente distribuidas = ambas siguen la misma distribución poblacional.
  • variables aleatorias IID = modelo por defecto para muestras aleatorias = punto inicial por defecto para el estudio de inferencia

VALOR ESPERADO O ESPERANZA MATEMATICA(Centro geométrico de las Áreas de probabilidad) o MEDIA

Se insiste mucho en la diferencia entre muestra y población. En este caso cuando hagamos cálculos sobre una muestra haremos los cálculos en base a los valores encontrados en esa muesta (Ej:tiro un dado 10 veces y salen 2 1s, 3 2s, 1 3, 4 5s.). Cuando se habla sobre los cáculos sobre una población estaremos hablando sobre las probabilidades que tenga una variable de tomar un valor sin haber realizado muestras ( o lo que es lo mismo si hubiéramos hecho infinitas muestras)

MEDIA DE UNA MUESTRA

Es la suma de todos los valores dividido por el número de muestras TirarDado = {1,1,5,4,3,4,6,5,4,6} _ x = (1+1+5+4+3+4+6+5+4+6)/10 = 3.9

MEDIA DE UNA POBLACIÓN

Es una propiedad de las distribuciones de probabilidad de la misma manera que lo son la mediana (valor del punto que está en la mitad) y la varianza

En variables discretas será el valor de su probabilidad multiplicado por el propio valor

E(x)= x1P(x1)+x2P(x2)+…

\[E[X] = 1 \times \frac{1}{6} + 2 \times \frac{1}{6} + 3 \times \frac{1}{6} + 4 \times \frac{1}{6} + 5 \times \frac{1}{6} + 6 \times \frac{1}{6} = 3.5\]

qownnotes-media-S12388

qownnotes-media-S12388

En variables contínuas, tenemos que hacer uso de la función de distribución para recoger la probabilidad acumulada.(dx)

E(x) = Integral (xf(x) d(x))

En PDF es lo mismo , el valor esperado es el centro de masa del área que encierra su función de densidad. Una propieda de la media es que su valor es el mismo en la función de densidad (DE POBLACION) que en la función de densidad de la media de X muestras del mismo tipo (MUESTRA) , es decir, si el valor esperado de una distribución normal centrada en 0 es 0 (en azul), función de densidad de la media de 10 distribuciones normales distintas será también 0 (en rosa). Cuanto mayor sea la población sobre la que estamos sacando la media (en lugar de 10 que fuera 100) más se aglutinarán las densidades entorno al valor esperado.

qownnotes-media-z14800

qownnotes-media-z14800

MEDIA DE MÚLTIPLES MUESTRAS

La media de un conjunto de muestras es en sí también una variable aleatoria, de la que podríamos calcular su media y su varianza ( la media de las medias y la varianza de las medias).

E(x) = suma x / n

VARIABILIDAD

VARIANZA DE UNA MUESTRA

Es la suma de la distancia de todos los puntos de la muestra respecto a la media (se pone al cuadrado para evitar distancias negativas y calcular su distancia absoluta)

S^2 = (1-3.9 +1-3.9 +5-3.9 +4-3.9 +3-3.9 +4-3.9 +6-3.9 +5-3.9 +4-3.9 +6-3.9 )^2 / 10 -1 = 0

S = Desviación estándar Se aplica a raiz cuadrada para que esté en las mismas unidades que el resto.

La varianza de una muestra se define como \[S^2 = \frac{\sum_{i=1} (X_i - \bar X)^2}{n-1}\]

También es una variable aleatoria:

  • Tiene asociada una distribución poblacional
  • Su valor experado es la varianza poblacional
  • Su distribución se concentra cada vez más alrededor de la varianza poblacional cuantos más datos haya.
  • Su raíz cuadrada es la desviación estándar

En el caso de la varianza de las muestras hay una pequeña variación ya que se divide por n-1. Esto es debido a que las primeras n-1 observaciones son son independientes de la media. la úiltima no es independiente ya que puede ser calculada por medio de la media de la muestra usada en la formula. En otras palabras, la varianza de la muestra es CASI la media de la desviación al cuadrado de la media de la muestra.

Var(x) suma var(x) /n-1

La distribución de la varianza de las muestras también se verá afectado por el tamaño de la muestra, de manera que a mayor tamaño de la muestra la curva se irá acercando a la variaza de la población que está intentando estimar:

qownnotes-media-M19368

qownnotes-media-M19368

Var(X’)=Var(1/n*Sum(X_i))=(1/n2)Var(Sum(X_i))=(1/n^2)Sum(sigma2)=sigma^2/n para poblaciones infinitas cuando las muestras son independientes | infinite populations when our samples are independent

VARIANZA DE UNA POBLACIóN

Varianza: es la dispersión de los datos respecto a la media.

Var(X) = E(X-mu)^2 (Mide la distancia de todos los puntos con respecto a la media, al cuadrado para que sean distancias absolutas) sd(X) = sqrt(var(x)) Desviación estándar es la raíz cuadrada de la varianza de manera que el dato está en las mismas unidades que el eje de las X, lo cual es muy útil.

Desviación típica = desviación estándar = error estándar = error típico

También tiene una distribución de población asociada así como un valor esperado llamado varianza de la población. La desviación estandar se llama también error estandar.

En el caso del ejemplo:Var(X) = E(x)^2 - mu^2

E(x)^2 = (1x1/6)2+(2x1/6)2+(3x1/6)2+(4x1/6)2+(5x1/6)2+(6x1/6)2 = 15.17

mu^2 = 3.5^2 = 12,25

Var(X) = 15.17-12.25 = 2.92

VARIANZA DE MÚLTIPLES MUESTRAS

La desviación estándar de una una estadística se llama error estándar y es la varianza dividido el número de muestras, o lo que es lo mismo la desviación estándar dividido raíz cuadrada del número de muestas (numerador y denominador aplicados la raíz cuadrada)

varianza = S^2/n desviación estándar = S/raizCuadrada(n)

MEDIA Y DESVIACIÓN ESTANDAR (ERROR ESTANDAR) DE MÚLTIPLES MUESTRAS EN DISTINTAS DISTRIBUCIONES

Si simulamos con R 1000 muestras de una distribción normal con 10 valores posibles vemos que se acerca mucho a lo calculado para la población, y cuanto más se aumenten las muestras más cerca estará ese valor:

Distribución Normal Estándar

# specify number of simulations with 10 as number of observations per sample
nosim <- 1000; n <-10
# estimated standard deviation of mean
sd(apply(matrix(rnorm(nosim * n), nosim), 1, mean))
## [1] 0.3037642
# actual standard deviation of mean of standard normals
1 / sqrt(n)
## [1] 0.3162278

Distribución estándar uniforme

  • Estandar uniforme \(\rightarrow\) triangle straight line distribution \(\rightarrow\) media = 1/2 and varianza = 1/12
  • Media de muestras aleatorias de \(n\) uniformes tienen desviación estándar de \(1/\sqrt{12 \times n}\)
# estimated standard deviation of the sample means
sd(apply(matrix(runif(nosim * n), nosim), 1, mean))
## [1] 0.09341847
# actual standard deviation of the means
1/sqrt(12*n)
## [1] 0.09128709

Distribucción de Poisson

  • \(Poisson(x^2)\) tiene varianza de \(x^2\)
  • Medias de muestras aleatorias de \(n~ Poisson(4)\) tienen una desviación estándar de \(2/\sqrt{n}\)
# estimated standard deviation of the sample means
sd(apply(matrix(rpois(nosim * n, lambda=4), nosim), 1, mean))
## [1] 0.6606151
# actual standard deviation of the means
2/sqrt(n)
## [1] 0.6324555

Distribución Binomial

  • Para \(p = 0.5\), La distribución de Bernuilli tiene varianza de 0.25
  • Medias de muestras aleatorias de \(n\) de valores 0 o 1 tienen una desviación estándar de \(1 / (2 \sqrt{n})\)
# estimated standard deviation of the sample means
sd(apply(matrix(sample(0 : 1, nosim * n, replace = TRUE), nosim), 1, mean))
## [1] 0.1601838
# actual standard deviation of the means
1/(2*sqrt(n))
## [1] 0.1581139

Hay una medida adicional que es el coeficiente de dispersión de Pearson, que es la desviación típica dividida por su media, de esta manera podemos medir la dispersión de dos variables de distintas unidades. Por ejemplo, comparar la dispersión del peso de los niños de una clase frente a la dispersión de la altura de otra clase.

El índice de Gini también nos puede servir para medir cómo de concentrados están los datos

El coeficiente de Fisher nos dice cómo de asimétrica es la curva de distribución.

La Courtosis es cómo de aglutinados están las frecuencias alrededor de la media.

EJEMPLOS CON DISTRIBUCIONES

BERNULLI

Si X es una variable aleatoria que mide el número de éxitos, y se realiza un único experimento con dos posibles resultados (éxito o fracaso), se dice que la variable aleatoria X se distribuye como una Bernoulli de parámetro p. (Ejemplo tirar 1 moneda 1 vez y que salga cruz)

  • Distribución de Bernoulli es el resultado de un experimento con sólo 2 resultados (binario) * 1 = “éxito” con probabilidad de \(p\) * 0 = “fallo” con probabilidad de \(1 - p\)
    • PMF se define como \[P(X=x) = p^x(1 - p)^{1-x}\]
    • mean = \(p\)
    • variance = \(p(1 - p)\)

Si hacemos varias observaciones de Bernuilli pongamos \(x_1,\ldots, x_n\) , la posibilidad (likelihood)de que el evento sea éxito es \[ \prod_{i=1}^n p^{x_i} (1 - p)^{1 - x_i} = p^{\sum x_i} (1 - p)^{n - \sum x_i} \] \(\hat p = \sum_i x_i / n\) es el máximo estimador de posibilidad para \(p\)

Todas las posibilidades para n pequeña

BINOMIAL

cuenta el número de éxitos en una secuencia de n ensayos de Bernoulli independientes entre sí, con una probabilidad fija p de ocurrencia del éxito entre los ensayos. Un experimento de Bernoulli se caracteriza por ser dicotómico, esto es, solo dos resultados son posibles

Para n = 1, la binomial se convierte, de hecho, en una distribución de Bernoulli.

  • variable aleatoria binomial = suma de n variables de Bernoulli \[X = \sum_{i=1}^n X_i\] donde \(X_1,\ldots,X_n = Bernoulli(p)\)
    • PMF se define como \[P(X=x) = {n \choose x}p^x(1-p)^{n-x}\] donde \({n \choose x}\) = número de maneras de seleccionar \(x\) veces de \(n\) opciones sin reemplazamiento \(x=0,\ldots,n\)
    • combinación o “\(n\) combinaciones de \(x\)” se define como \[{n \choose x} = \frac{n!}{x!(n-x)!}\]
    • Los casos base son \[{n \choose n} = {n \choose 0} = 1\]

qownnotes-media-dnOuDU FuncionMasaProbabilidad

DistribucionAcumulada * de 8 hijos, cuál es la probabilidad de tener 7 chicas o más (50/50 chance)? \[{8 \choose 7}.5^7(1-.5)^{1} + {8 \choose 8}.5^8(1-.5)^{0} \approx 0.04\]

  • choose(8, 7) = función de R function para calcular las combinaciones de x en n elementos
  • pbinom(6, size=8, prob =0.5, lower.tail=TRUE) = probabilidad de 6 o menos éxitos de 8 muestras con probabilidad 0.5 (CMF)
    • lower.tail=FALSE = devuelve el complemento, en este caso es la probabilidad de que sea mayor a 6 de 8 muestrasy of 0.5
qownnotes-media-fjsoAv

qownnotes-media-fjsoAv

plot(pvals, dbinom(7, 8, pvals) / dbinom(7, 8, 7/8) ,
     lwd = 3, frame = FALSE, type = "l", xlab = "p", ylab = "likelihood")

DISTRIBUCION NORMAL ESTANDAR

Es una distribución normal pero con la particularidad de que tiene mu = 0 y sigma = 1

qownnotes-media-BjqEKS

qownnotes-media-BjqEKS

De manera que 1 desviación estandar recoge el 68% de los datos, 2 desviaciones estándar el 95% y 3 desviaciones estandar el 99%.

qownnotes-media-SIFoSu

qownnotes-media-SIFoSu

qownnotes-media-hjbGgE

qownnotes-media-hjbGgE

# specify number of simulations with 10 as number of observations per sample
nosim <- 1000; n <-10
# estimated standard deviation of mean
sd(apply(matrix(rnorm(nosim * n), nosim), 1, mean))
## [1] 0.3166805
# actual standard deviation of mean of standard normals
1 / sqrt(n)
## [1] 0.3162278

DISTRIBUCION NORMAL

Por tanto una distribucion normal tendrá mu != 0 y sigma != 1, la única diferencia con los gráficos anteriores es que no estará centrada en 0, y que los ejes de las x cada 1,2 y 3 estar? multiplicado por su sigma.

DistribuciónNormal

DistribuciónNormal

qownnotes-media-mAfrpi

qownnotes-media-mAfrpi

  • Approximadamente en 68%,95% y 99% de la densidad normal recae en 1,2 y 3 desviaciones estándar de la media.

  • -1,28,-1.645,-1.96 and -2.33 son los percentiles 10º,5º, 2.5º y 1º de la distribución

  • 1,28,1.645,1.96 and 2.33 son por simetría los percentiles 90º,95º, 97.5º y 99º de la distribución

qownnotes-media-RxbLpG

qownnotes-media-RxbLpG

qownnotes-media-R14800

qownnotes-media-R14800

qownnotes-media-B14800

qownnotes-media-B14800

qownnotes-media-P19368

qownnotes-media-P19368

qownnotes-media-L19368

qownnotes-media-L19368

POISSON

expresa, a partir de una frecuencia de ocurrencia media, la probabilidad de que ocurra un determinado número de eventos durante cierto periodo de tiempo. Concretamente, se especializa en la probabilidad de ocurrencia de sucesos con probabilidades muy pequeñas, o sucesos “raros”.

Tanto el valor esperado como la varianza de una variable aleatoria con distribución de Poisson son iguales

qownnotes-media-mMWyDG

qownnotes-media-mMWyDG

qownnotes-media-BmoCZj

qownnotes-media-BmoCZj

qownnotes-media-Y14800

qownnotes-media-Y14800

qownnotes-media-i14800

qownnotes-media-i14800

qownnotes-media-i19368

qownnotes-media-i19368

qownnotes-media-u19368

qownnotes-media-u19368

Es muy útil para:

*Conteos de datos Cuantas veces se repite un determinado evento.

*Modelar eventos de tiempo o de supervivencia, que modeliza el tiempo que se tarda en que ocurra un determinado suceso. Por el nombre de la técnica parecería que se analizara el tiempo hasta la muerte (Análisis de supervivencia) pero, en realidad, puede analizarse cualquier otro suceso.

En Análisis de supervivencia la muestra consiste en el seguimiento de una serie de individuos desde el inicio del estudio hasta su final y, ante una situación de este tipo, es frecuente que se produzca la desaparición de alguno de esos individuos que entran en el estudio. También es posible que al entrar un individuo en el estudio, éste termine antes de que en ese individuo se produzca el suceso que se pretende detectar. Un ejemplo es el tiempo en que cominezan a aparecer determinados sintomas en pacientes sanos: qownnotes-media-gtvKdN

*Para modelar tablas de contingencia.Modelando tablas de contingencia (cruzar dos características y ver el número de casos que pertenecen a cada categoria)

qownnotes-media-uSdaBj

qownnotes-media-uSdaBj

*Aproximaciones binomiales cuando n es muy grande y p muy pequeño, por ejemplo en epidemiología el tener o no tener una enfermedad sobre una población.

qownnotes-media-p14800

qownnotes-media-p14800

qownnotes-media-l14800

qownnotes-media-l14800

qownnotes-media-n19368

qownnotes-media-n19368

qownnotes-media-m19368

qownnotes-media-m19368

ESTUDIO ASINTÓTICO

Es el término para el análisis del comportamiento de las estadísticas cuando el tamaño de la muestra tiene un límite (normalmente infinito). Es la base de estudio frecuentista de probabilidades.

Un estimador es consistente si converge con lo que queremos estimar. La ley de números largos (LLN) dice que la media de la muestra, la varianza y la desviación estandar de iid muestras es consistente con la media de la población, la varianza y la desviación estándar respectivamente.

Teorema central del limite (CLT) : dice que la distribución de las medias siguen una distribución normal estándar centrada en el origen (ya que estamos restando la media de la población), con una desviación estándar igual al error estándar de la media, pero debe ser para un número de repeticiones suficiente como para que se cumpla.

\(X_1,\ldots,X_n\) es una colección de variables aleatorias iid con media \(\mu\) y varianza \(\sigma^2\) \(\bar X_n\) es la media de la muestra Entonces \(\frac{\bar X_n - \mu}{\sigma / \sqrt{n}}\) tiene una distribución normal para largos \(n\).

Por ejemplo, si tiramos un dado muchas veces y hacemos n observaciones

\(\frac{\bar X_n - 3.5}{1.71/\sqrt{n}}\)

Se ve claramente la dierencia de hacer 1,2 y 6 tiradas 10000 veces. La curva normal estándar lógicamente va a estar centrada en el origen ya que le estamos restando la media.

par(mfrow = c(1, 3))
for (n in c(1, 2, 6)){
  temp <- matrix(sample(1 : 6, n * 10000, replace = TRUE), ncol = n)
  temp <- apply(temp, 1, mean)
  temp <- (temp - 3.5) / (1.71 / sqrt(n))
  dty <- density(temp)
  plot(dty$x, dty$y, xlab = "", ylab = "density", type = "n", xlim = c(-3, 3), ylim = c(0, .5))
  title(paste("Media de la muestra con ", n, "observaciones"))
  lines(seq(-3, 3, length = 100), dnorm(seq(-3, 3, length = 100)), col = grey(.8), lwd = 3)
  lines(dty$x, dty$y, lwd = 2)
}

Al seguir una distribución normal, podemos usar los cuantiles para saber con qué probabilidad la media obtenida estará incluida en ese intervalo, para una distrióbucin normal sabemos \[\bar X_n \pm z_{1-\alpha/2}\sigma / \sqrt{n}\] y lo llamaremos intervalo de confianza.

En este ejemplo de las alturas de los hijos, podemos dar un intervalo de confianza en el cuál con un 95% de posibilidades la media obtenida en la muestra está contenida en la media de la población aplicando esta ecuación

library(UsingR);data(father.son);
x <- father.son$sheight
#Intervalo de confianza calculado a mano
(mean(x) + c(-1, 1) * qnorm(.975) * sd(x) / sqrt(length(x))) / 12
## [1] 5.709670 5.737674

Esta misma operación la podemos hacer con Bernuilli para proporciones de éxito o fallo:

\[ \hat p \pm z_{1 - \alpha/2} \sqrt{\frac{p(1 - p)}{n}} \] Sabiendo que la moneda es imparcial, p = 1/2, p(1-p)=1/4 y sqrt(p(1-p)) = 1/2

Aplicando le intervalor Wald del 95%

\(p(1-p) \leq 1/4\) for \(0 \leq p \leq 1\) \(\alpha = .05\) so that \(z_{1 -\alpha/2} = 1.96 \approx 2\) \[ 2 \sqrt{\frac{p(1 - p)}{n}} \leq 2 \sqrt{\frac{1}{4n}} = \frac{1}{\sqrt{n}} \] Entonces \(\hat p \pm \frac{1}{\sqrt{n}}\) es un buen intervalo de confianza rápido para \(p\)

Un ejemplo, es que en un sondeo de 100 personas, 56 votan sí, tengo la votación ganada? 1/sqrt(100)=.1 esto nos da un intervalo de confianza de Wald (95%) de entre (0.46, 0.66) así que no tenemos nada claro.

(0.56 + c(-1, 1) * qnorm(.975) * sqrt(0.56*0.44/100))
## [1] 0.4627099 0.6572901
binom.test(56,100)$conf.int
## [1] 0.4571875 0.6591640
## attr(,"conf.level")
## [1] 0.95

Si el proceso no fuera parcial (p=0.9 por ejemplo) el CTL nos garantiza únicamente que es una buena aproximación en el límite, pero se puede observar que si hacemos la misma prueba que con p = 0.5 tardan muchas más iteraciones en converger hacia una normal estándar.

Hay veces que el intervalo de wald no se comporta de manera correcta si el número de repeticiones no es el adecuado, y para evitarlo se puede aplicar la corrección de Agresti/Coull añadiendo dos aciertos y dos errores $p = X_n +2 / n+4 $

Por último si se trata de una distribución de Poisson por ejemplo para conteos:

Vemos que damos un un error diagnóstico 5 veces cada 94.32 días, cuál es el porcentaje de fallo diario con un intervalo de confianza dentro del 95%?

\(X \sim Poisson(\lambda t)\) Estimate \(\hat \lambda = X/t\) \(Var(\hat \lambda) = \lambda / t\)

\[ \frac{\hat \lambda - \lambda}{\sqrt{\hat \lambda / t}} = \frac{X - t \lambda}{\sqrt{X}} \rightarrow N(0,1) \]

x <- 5
t <- 94.32
lambda <- x/t
round(lambda + c(-1, 1) * qnorm(.975) * sqrt(lambda/t),3 )
## [1] 0.007 0.099
poisson.test(5,94.32)$conf.int
## [1] 0.01721254 0.12371005
## attr(,"conf.level")
## [1] 0.95

En este caso podemos observar que el intervalo calculado es más conservador ya que t no es lo suficientemente grande.

qownnotes-media-bG7876

qownnotes-media-bG7876

Hay veces que el número de repeticiones no es suficiente en distribuciones discretas binomiales y de Poisson y se puede hacer una aproximación del valor estimado sumando +2 al numerador y +4 al denominador (añadiendo 2 exitos y 2 fallos)

qownnotes-media-jY7876

qownnotes-media-jY7876

qownnotes-media-ZUuUAB

qownnotes-media-ZUuUAB

INTERVALOS DE CONFIANZA

Resumiendolo que hemos visto anteriormente, la media de los valores medios se aproxima a una normal estandar, con media \(\mu\) y sd \({\sigma / \sqrt{n}}\).

Por tanto, al seguir una normal estandar, el desplazamiento desde la media de 2 unidades del error estándar nos dará un intervalo en el que tendremos la probabilidad del 95% de contener el valor medio. A este intervalo lo llamamos ZQuantile X SE

\(\bar X_n \pm 2\sigma / \sqrt{n}\) es el llamado intervalo 95% de \(\mu\)

En la vida real, las muestras con las que trabajamos no suelen ser tan grandes, por lo que necesitamos una adaptación, utilizando la distribución T Student de Gosset. surge del problema de estimar la media de una población normalmente distribuida cuando el tamaño de la muestra es pequeño. Está también centrada en 0. Si no fuera simétrica en 0 (skewed distr) habría que calculalo de otra manera o coger otro estimador como la mediana

qownnotes-media-WYnisx

qownnotes-media-WYnisx

Aparece de manera natural al realizar la prueba t de Student para la determinación de las diferencias entre dos medias muestrales y para la construcción del intervalo de confianza para la diferencia entre las medias de dos poblaciones cuando se desconoce la desviación típica de una población y ésta debe ser estimada a partir de los datos de una muestra.

qownnotes-media-MIqXgZ

qownnotes-media-MIqXgZ

En adelante vamos a ver pequeñas muestras de datos y en este caso tendremos que usar la distribución T de Gosset y intervalos de confianza T, la única diferencia con la ecuación anterior es que usaremos el TQuantile

Cuando una distribución T creece y crece en muestras se acaba convirtiendo en distribución Z.

La distribución T tiene los extremos más pequños que la normal.A diferencia de la dist estandar, que tiene 2 parametros (media y varianza). La T de gosset tiene sólo un parámetro llamado grados de libertad t = n-1, cuanto más crezca este número más se parecerá a una standar normal, ya que si n es muy grande la diferencia entre una estandar normal, que es dividida por sigma, y la T de student que es dividida por la varianza, se hace insignificativa, es algo similar a lo que nos pasaba al estudiar la varianza de las muestras, que cuanto mayor era el número de las muestras más se acerca al teórico de la población. En este caso para calcular T con n-1 grados de libertad:

qownnotes-media-KN1524

qownnotes-media-KN1524

El valor de esta t nos dará un intervalo de confianza, para el cuál dependiendo del número de individuos, en 95% de los casos ese valor está incluido, en el caso de que quisiéramos usar ese quantil qt(0.975,n)

Dependiendo de los grados de confianza, la diferencia de las distribuciones será mayor, en este caso los grados de confianza son muy amplios, lo que hace que sea muy similar

qownnotes-media-ASAszo

qownnotes-media-ASAszo

La distribución T es muy útil cuando queremos comparar cómo varía un grupo frente a otro, pero siempre teniendo en cuenta que deben ser datos iid normales, centrados en su media

Vamos a coger dos grupos y vamos a estudiar la diferencia de variabilidad entre ellos, primero estando relaccionados (que los individuos sean los mismos, y queremos ver su evolución) y luego siendo grupos independientes (Son grupos de personas totalmente distintas) .

Esa diferencia de variabilidad estará encerrada en un intervalo de confianza y se suele utilizar el intervalo T para compararlos

Los resultados del estudio por tanto varían dependiendo de si se emparejan los resultados porque están relaccionados o son independientes: Los datos son un grupo de 10 pacientes a los que se les ha tratado con 2 drogas distintas y se intenta comparar los resultados

data(sleep)
sleep
##    extra group ID
## 1    0.7     1  1
## 2   -1.6     1  2
## 3   -0.2     1  3
## 4   -1.2     1  4
## 5   -0.1     1  5
## 6    3.4     1  6
## 7    3.7     1  7
## 8    0.8     1  8
## 9    0.0     1  9
## 10   2.0     1 10
## 11   1.9     2  1
## 12   0.8     2  2
## 13   1.1     2  3
## 14   0.1     2  4
## 15  -0.1     2  5
## 16   4.4     2  6
## 17   5.5     2  7
## 18   1.6     2  8
## 19   4.6     2  9
## 20   3.4     2 10

En las siguientes gráficas se ven cada uno de los valores de ambos resultados, estando la leyenda la correspondencia entre el color de la raya que une a los individuos y el número de individuo que es en el grupo

DEPENDIENTES (paired=TRUE)

## [1] -3.363874  0.203874

qownnotes-media-JY1524

qownnotes-media-JY1524

qownnotes-media-rV1524

qownnotes-media-rV1524

Si queremos calcular el intevalo de confianza de los resultados emparejados, creamos un vector nuevo con la diferencia de los resultados de ambos grupos difference <- g2-g1

y se calcula la media y la desviación estandar mn <- mean(difference) s <- sd(difference) Los calculos se pueden hacer a mano mn + c(1,-1)qt(.975)(s/sqrt(10)) o con la función t.test()

qownnotes-media-qs1524

qownnotes-media-qs1524

qownnotes-media-ay1524

qownnotes-media-ay1524

Esto lo que nos dice, es que con 95% de posibilidades ( la funcion qt(.975) nos da la t para el 95% con 9 grados de libertad) la diferencia media de los efectos entre ambas drogas es entre 0.7 y 2.46 horas de sueño

data(sleep)
x1 <- sleep$extra[sleep$group == 1]
x2 <- sleep$extra[sleep$group == 2]
n1 <- length(x1)
n2 <- length(x2)
sp <- sqrt( ((n1 - 1) * sd(x1)^2 + (n2-1) * sd(x2)^2) / (n1 + n2-2))
md <- mean(x1) - mean(x2)
semd <- sp * sqrt(1 / n1 + 1/n2)
md + c(-1, 1) * qt(.975, n1 + n2 - 2) * semd
## [1] -3.363874  0.203874
t.test(x1, x2, paired = FALSE, var.equal = TRUE)$conf
## [1] -3.363874  0.203874
## attr(,"conf.level")
## [1] 0.95
t.test(x1, x2, paired = TRUE)$conf
## [1] -2.4598858 -0.7001142
## attr(,"conf.level")
## [1] 0.95

INDEPENDIENTES (PAIRED= FALSE)

En este caso tenemos dos grupos distintos de personas, unos que tomaron el medicamento y otro que tomaron placebo y vamos a medir su presión arterial. No se pueden emparejar porque son personas distintas y los grupos pueden tener incluso tamaño diferente.

qownnotes-media-Fz1524

qownnotes-media-Fz1524

El objetivo de nuevo es encontrar el intervalo de confianza de la diferencia media de la presión arterial media entre ambos grupos con un 95% de posibilidad. Para hacer el cálculo de los intervalos de confianza vamos a trabajar restando la media de ambos conjuntos, el quantile se calcula en base a ambos conjuntos, ya que hemos calculado su media por separado, así que hay que restar 2 a la suma de ambos para obtener el grado de confianza de todo.

qownnotes-media-aV1524

qownnotes-media-aV1524

La varianza común es más dificil de calcular, haremos una estimación media entre ambas varianzas aplicando pesos en cada una de ellas dependiento del tamaño de cada grupo diviendolo por el tamaño total. Se ve más claro en este ejemplo extraído del libro Fundamentals of Biostatistics. Primero haremos los cálculos suponiendo que las varianzas son equivalentes. Se puede observar que el intervalo de confianza sobre difrencia media de ambas presiones arteriales medias contiene el valor 0, por lo que no se puede descartar que las medias de ambos grupos sean iguales.

qownnotes-media-Tn1524

qownnotes-media-Tn1524

Aquí se puede apreciar el error de cálculo si tratamos el ejemplo anterior como grupos independientes

qownnotes-media-kP1524

qownnotes-media-kP1524

Se puede apreciar en el cálculo de ambos casos que cuando los tratamos de manera distinta un cálculo contiene el 0 y el otro no lo contiene, lo que significa que cuando los grupos están emparejados tienen muchas más medias distintas, de manera opuesta a cuando no los emparejamos y los tratamos como independientes en que se contiene el valor 0 y por tanto las diferencias en las medias son menos distintas ya que se tratan todos los datos como un solo grupo.

La función t.test tiene otro parámetro con el que se puede jugar que es si la varianza de los datos es igual en ambos conjuntos. En este ejemplo, se intenta hacer un intervalo de confianza entre la ganancia de peso de pollos con 4 dietas, comparándolas de 2 en 2, en este ejmplo la 1 y la 4, cuya diferencia de variablidad es muy evidente a nivel gráfico.

qownnotes-media-Bt1524

qownnotes-media-Bt1524

qownnotes-media-mA1524

qownnotes-media-mA1524

vemos que los resultados de las estimaciones varían dependiendo del tratamiento de la relacción entre las varianzas de los grupos.

qownnotes-media-ll1524

qownnotes-media-ll1524

DISTINTA VARIANZA ENTRE GRUPOS TEST DE WELCH

En el caso de que la varianza de los grupos no sea la misma o no estemos seguros, En caso de duda, se debe utilizar varianza no es igual, para calcularala

qownnotes-media-Mp1524

qownnotes-media-Mp1524

qownnotes-media-Rh1524

qownnotes-media-Rh1524

CHI SQUARE TEST (DATOS CATEGÓRICOS)

Este test se utiliza para saber si dos variables son independientes, es decir, si la aparición de una está relaccionada (ojo no es la causa, sólo si están relacconada) con la otra. Es apropiado utilizar este método en las siguientes situaciones:

  • El muestreo es aleatorio
  • las variables son categoricas
  • si los datos recogidos en la tabla de contingencia tienen en las frecuencias esperadas al menos 5 casos.

Este método consiste en 4 pasos:

Establecer la hipótesis: Ho-> Variable A y variable B son independientes Ha-> Variable A y variable B no son independientes . La hipótesis alternativa por tanto es que sabiendo el valor A nos puede ayudar a predecir el valor de B Formular el plan de análisis Este paso describe cómo usar muestras de los datos para aceptar o rechazar la hipótesis nula: Nivel de significancia elegir un nivel apropiado entre 0 y 1 Test a utilizar en este caso chi-cuadrado.

Analizar la muestra de los datos Grados de libertad DF = (r-1)(c-1) donde r son el nº de niveles de A y c el nº de niveles de B Frecuencia esperada = Son computadas por separado para cada nivel de cada variable en cada nivel de la otra variable Er,c = (nr*nc)/n Donde Er,c es la frecuencia esperada para el nivel r de la variable A y nivel c de la B nr es el número total de observaciones del nivel r en A nc es el número total de observaciones del nivel c en B n es el número total de la muestra Test estadístico = x2 = sum[(Or,c - Er,c)2/Er,c] Donde Or,c es la frecuencia del nivel r de la variable A y nivel C de la B Er,c es la frecuencia esperada del nivel r de la variable A y el nivel c de la B P-value = es la probabilidad de observar una muestra de manera aleatoria al menos tan extrema como la nuestra. Para poder conocer esta probabilidad, se debe calcular si la distribución pertenece a una chi-cuadrado, si es así, utilizar los grados de libertad anteriores y el test estadístico para obtener el pvalue de la muestra de datos.

CONTRASTE DE HIPÓTESIS (Hypothesis testing) DATOS CON DISTRIBUCIONES NORMALES

Un contraste de hipótesis (también denominado test de hipótesis o prueba de significación) es un procedimiento para juzgar si una propiedad que se supone en una población estadística es compatible con lo observado en una muestra de dicha población.La aplicación de cálculos probabilísticos permite determinar a partir de qué valor debemos rechazar la hipótesis garantizando que la probabilidad de cometer un error es un valor conocido a priori.

En el contraste de hipótesis tenemos sempre dos resultados:

Una hipotesis nula: H0 (de no efecto) Una hipotesis alternativa: H1 (de diferencia o efecto)

A partir de una muestra de la población en estudio, se extrae un estadístico (media muestral, varianza muestral,etc.) cuya distribución de probabilidad esté relacionada con la hipótesis en estudio y sea conocida. Se toma entonces como región de rechazo al conjunto de valores que es más improbable bajo la hipótesis, esto es, el conjunto de valores para el que rechazaremos la hipótesis nula si el valor del estadístico observado entra dentro de él.

Para probar una hipótesis vamos a comparar siempre con la hipótesis nula Ho, e invalidarla, ya que partimos de que la hipótesis nula es siempre cierta.

A continuación, vamos a ver cuatro ejemplos de cómo descartar H_0, con una distribución normal grande, con una distribución normal pequeña ambos con varianza equvalente, uno con varianza no equivalente y un ejemplo con distribición binomial.

Z-TEST

Vamos a trabajar con un conjunto de 100 pacientes con sobrepeso que tienen desorden del sueño y una media de 32 eventos/hora durante la noche con una desviación estándar de 10 eventos/hora. Se considera que desorden del sueño es 30 eventos/hora. Vamos a utilizar el Test-Z como test estadístico, para poder usarlo necesitamos la media , la varianza y la distribución.

Queremos invalidar H0 que dice que mu = 30 y validar H_a que dice que mu > 30 (32 en este caso concreto)

ERROR TIPO I –> EL investigador SI encuentra diferencia entre Ho y Ha cuando no la hay (Falso positivo) ALFA probabilidad de cometer ese error

ERROR TIPO II –> El investigador NO encuentra diferencia entre Ho y Ha cuando sí la hay (falso negativo) BETA probabilidad de cometer ese error

Error tipo I –> Rechazar Ho cuando es correcta esto es alfa (probabilidad de cometer este error es típicamente 5% llamada alfa) Tampoco queremos un alfa muy bajo ya que nunca rechazaríamos H_o aunque fuera falsa. Error tipo II –> Aceptar H_o cuando es falso Error tipo III –> Rechazar H_a cuando es correcta Error tipo IV –> Aceptar H_a cuando es falsa

qownnotes-media-amRzlL

qownnotes-media-amRzlL

El error estándar de la media de un conjunto es desviacion_estandar/sqrt(población), en este ejemplo 10/sqrt(100) = 1

H_o por tanto tiene un media normalmente distribuida,mu = 30 y varianza = 1 (estimada como el error estandar que es 1 al cuadrado)

Estamos buscando una constante C que hiciera que la probabilidad de que la media sea mayor que C es %5 Esto es P(X>C|H_o) = 5%

!(media/32517.png)

El área sin colorear se puede calcular como qnorm(.95)

Recordando que en una distribución normal el percentil 95 está a 1.645 desviaciones estándar de la media y teniendo una distribución normal de media 30 y 1 de desviación estandar N(30,1)

Rechazaríamos la hipotesis nula si la media fuera mayor que un valor (C) cuya probabilidad sea 5% C = 30 + 1x1.645 = 31.645

!(media/26308.png)

Rechazar la hipótesis nula cuando la media es >= 31.645 tiene una probabilidad de fallo del 5%, nuestra media es 32 que es mayor y cae dentro de ese área así que rechazamos H_o En este gráfico vemos las dos distribuciones, y cómo la media de H_a está dentro del área del 5% de H_o, lo que nos sirve para rechazar que nuestro resultado sea por azar en un 95% de los casos al menos.

!(media/26001.png)

El Z-score es lo mismo pero calculando cuantos errores estándar la media de la hipótesis es diferente de la media, esto se calcula computando la distancia entre las medias (32-30) y dividiendo por el error estandar de la media (s/sqrt(n)) que en este caso es 10/sqrt(100)

Z-score = 2, significa que está a 2 errores estandar de la media.

Un alumno A saca una puntuación de 85 en un examen cuyas puntuaciones tienen una media de 79 con una desviación típica de 8. Un alumno B saca 74 en un examen cuyas puntuaciones tienen una media de 70 y desviación típica de 5. ¿Cuál de los dos alumnos obtuvo una puntuación mejor? La respuesta, desde el punto de vista de la unidad tipificada, se obtiene así:

Las puntuaciones tipificadas de los alumnos A y B son respectivamente:

qownnotes-media-blJJQy

qownnotes-media-blJJQy

Así el alumno B lo hizo mejor que el A, aunque su puntuación de 74 es inferior a 85.

qownnotes-media-rn7296

qownnotes-media-rn7296

La regla por tanto para descartar H_0 es (X-mu)*sqrt(n)/s > Z

T-TEST

Digamos que ahora bajamos de 100 ejemplos a 16, se convierte entonces en una distribución T con 15 grados de libertad. Aquí ya no podemos utlizar Zscore como test estadístico porque la distribución es normal, pero con pocos individuos (distribución T-student)

qownnotes-media-Hy7296

qownnotes-media-Hy7296

En este caso fallamos al rechazar Ho porque 0.8< 1.7531

En lugar de concentrar todo el error en uno de los extremos, se suele distribuir a ambos lados

qownnotes-media-sv7296

qownnotes-media-sv7296

Normalmente todos estos cálculos no los tenemos que hacer a mano, ya que la función t.test() nos da toda esta información

qownnotes-media-nS7296

qownnotes-media-nS7296

En este ejemplo, queremos saber si la media de altura de los hijos es equivalente a la media de altura de los padres. Cogemos los datos de las alturas de padres e hijos, y restamos la altura de un padre a la de su hijo, y vamos a ver si la diferencia (haciendo la resta se asume que están emparejados) es 0 o no es 0 con t.test (tambien se podrían pasar ambos arrays a la función y como tercer parametro paired.equals=TRUE)

Ambos t(cuantas desviaciones estandar está la media del conjunto del de H_o) y df (grados de libertad) son muy grandes con lo que no hace falta calcularlo.

Otro ejemplo

qownnotes-media-cy7296

qownnotes-media-cy7296

Effect Size

El tamaño del efecto o también llamado effect size es la diferencia neta observada entre los grupos de un estudio. No se trata de una medida de inferencia estadística ya que no se pretende identificar si las poblaciones son significativamente diferentes, sino que simplemente indica la diferencia observada entre muestras, independientemente de la varianza que tengan. Se trata de un parámetro que siempre debe acompañar a los p-values, ya que un p-value solo indica si hay evidencias significativas para rechazar la hipótesis nula pero no dice nada de si la diferencia es importante o práctica. Esto último se averigua mediante el tamaño del efecto.

En el caso de los t-test de medias independientes, existen dos medidas posibles del tamaño del efecto: la d de Cohen y la r de Pearson. Ambas son equivalentes y pueden transformarse de una a otra. Cada una de estas medidas tiene unas magnitudes recomendadas para considerar el tamaño del efecto como pequeño, mediano o grande. La función cohen.d() del paquete effsize permite calcular el tamaño del efecto de la diferencia de medias independientes.

D DE COHEN PARA T-TEST INDEPENDIENTES Y DEPENDIENTES: sd= sqr(sd2a+sd2b/na+nb) Los límites más utilizados para clasificar el tamaño del efecto con d-Cohen son:

d < 0.2 pequeño d > 0.5 mediano d = 0.8 grande

library(effsize)
datos <- data.frame(corredor = c(1:10),
                    antes = c(12.9, 13.5, 12.8, 15.6, 17.2, 19.2, 12.6, 15.3,
                              14.4, 11.3),
                    despues = c(12.7, 13.6, 12.0, 15.2, 16.8, 20.0, 12.0, 15.9,
                                16.0, 11.1))
cohen.d(d = datos$antes, f = datos$despues, paired = TRUE)
##
## Cohen's d
##
## d estimate: -0.06745406 (negligible)
## 95 percent confidence interval:
##     lower     upper
## -1.007282  0.872374

R DE PEARSON PARA T-TEST INDEPENDIENTE (NO VALE PARA DEPENDIENTES):

r=sqr(t^2*(t^2+gl))

t = estadístico t obtenido en el test gl = grados de libertad del test Los límites más utilizados para clasificar el tamaño del efecto con r son:

d < 0.1 pequeño d > 0.3 mediano d = 0.5 grande

BINOMIAL

En este caso binomial con valores discretos, vemos que es imposible coger el 5%, lo más cercano sería 0.0352, no podremos elegir ni Zscore ni Ttest, por lo que tendremos que hechar mano del PValue

qownnotes-media-iN7296 qownnotes-media-wd7296

CONTRASTE DE HIPÓTESIS (Hypothesis testing) NO NORMALES

A la hora de hacer contraste de hipótesis, podemos hablar de un contraste cuantitativo (cuán diferentes es una ditribución de algo medido en datos CONTINUOS, por ejemplo días de estancia, respecto de otra, ya sea su distribución supuesta u otra) o un contraste cualitativo (basado en el conteo de ocurrencias de variables categóricas DISCRETOS, por ejemplo tipos de altas al mes, si dos conteos tienen diferencias significativas, ya sea frente a su distribución supuesta u otra).

Otro factor fundamental a la hora de hacer el contraste será si conocemos su distribución supuesta (Test paramétricos, porque usaremos sus parámetros poblacionales media y varianza) o si no podemos definirla (No paramétricos).

Para hacer esta primera valoración, miraremos si cumple los criterios de normalidad y si sus varianzas son equivalentes.

Test de Normalidad Shapiro-Wilk Este test se emplea para contrastar normalidad cuando el tamaño de la muestra es menor de 50. Para muestras grandes es equivalente al test de kolmogorov-Smirnov. Uno de los supuestos que debemos de tratar en este tipo de prueba es el de distribucion normal. Utilizaremos para ello el estadistico de Shapiro-Wilk, existen otros estadisticos que pueden tambien utilizarse para probar la distribucion normal.El problema que nos encontramos es que para más de 5000 elementos este método ineficaz. Test de Kolmogorov-Smirnov Este es otro tipo de prueba que se utiliza para comparar la normalidad de un conjunto de datos. La prueba determina si los datos provienen de una misma distribucion. Puede ser utilizado tambien para comparar si un conjunto de datos se ajustan a una distribucion (en nuestro caso normal).

Cuando las muestras son pequenas la prueba de Kolmogorov es una prueba menos poderosa para detectar diferencias significativas.

El test de Kolmogorov-Smirnov permite estudiar si una muestra procede de una población con una determinada distribución (media y desviación típica), no está limitado únicamente a la distribución normal. Se ejecuta con la función ks.test().

ks.test(x = mtcars$mpg,"pnorm", mean(mtcars$mpg), sd(mtcars$mpg))
## Warning in ks.test(x = mtcars$mpg, "pnorm", mean(mtcars$mpg),
## sd(mtcars$mpg)): ties should not be present for the Kolmogorov-Smirnov test
##
##  One-sample Kolmogorov-Smirnov test
##
## data:  mtcars$mpg
## D = 0.1263, p-value = 0.687
## alternative hypothesis: two-sided

A pesar de que continuamente se alude al test Kolmogorov-Smirnov como un test válido para contrastar la normalidad, esto no es del todo cierto. El Kolmogorov-Smirnov asume que se conoce la media y varianza poblacional, lo que en la mayoría de los casos no es posible. Esto hace que el test sea muy conservador y poco potente. Para solventar este problema, se desarrolló una modificación del Kolmogorov-Smirnov conocida como test Lilliefors. El test Lilliefors asume que la media y varianza son desconocidas, estando especialmente desarrollado para contrastar la normalidad. Es la alternativa al test de Shapiro-Wilk cuando el número de observaciones es mayor de 50. La función lillie.test() del paquete nortest permite aplicarlo.

library("nortest")
lillie.test(x = mtcars$mpg)
##
##  Lilliefors (Kolmogorov-Smirnov) normality test
##
## data:  mtcars$mpg
## D = 0.1263, p-value = 0.2171

Respecto a la prueba de normalidad de Shapiro, esta resulta ser mas sensitiva que Kolmogorov. Es probable que la primera muestre que no existe normalidad y la segunda que si.

P1<-c(111,128,111,117,119,90,115,118,109,115,118,113,92,79,75,93,86,99,118)
P2<-c(88,118,100,107,119,87,105,121,111,117,108,103,96,82,72,90,91,92,115)
shapiro.test(P1)
##
##  Shapiro-Wilk normality test
##
## data:  P1
## W = 0.89668, p-value = 0.04239
mean<-mean(P1)
sd<-sd(P1)
ks.test(P1, "pnorm", mean=mean(P1), sd=sd(P1))
## Warning in ks.test(P1, "pnorm", mean = mean(P1), sd = sd(P1)): ties should
## not be present for the Kolmogorov-Smirnov test
##
##  One-sample Kolmogorov-Smirnov test
##
## data:  P1
## D = 0.21985, p-value = 0.3174
## alternative hypothesis: two-sided

Como observa, en este caso Shapiro muestra que no existe normalidad de datos (P1), mientras que Kolmogorov si lo muestra, en el mismo vector, demostrando asi la sensibilidad de cada una de las pruebas.

ks.test(P1,P2)
## Warning in ks.test(P1, P2): cannot compute exact p-value with ties
##
##  Two-sample Kolmogorov-Smirnov test
##
## data:  P1 and P2
## D = 0.31579, p-value = 0.2997
## alternative hypothesis: two-sided

Test de Jarque-Bera El test de Jarque-Bera no requiere estimaciones de los parámetros que caracterizan la normal. El estadístico de Jarque-Bera cuantifica que tanto se desvían los coeficientes de asimetría y curtosis de los esperados en una distribución normal.

Los datos presentan normalidad. Recuerde siempre visualizar de forma grafica los datos, tanto por su función de densidad (PDF) como por su función de distribución (CDF)

par(mfrow=c(2,2))
hist(P1)
hist(P2)
plot(ecdf(P1),xlim=range(c(P1,P2)),main="Distribucion Empirica Acumulada\n(P1,P2)")
plot(ecdf(P2), add=TRUE, lty="dashed", col="red")

Tipos de distribuciones

par(mfrow=c(2,2))
qqnorm(P1,col="red")
qqline(P1)
qqnorm(P2,col="red")
qqline(P2)
qqplot(P1,P2, main="Q-Q Plot de P1 y P2")

Si las gráficas tienen la misma distribución, el gárfico Q-Q entre grupos mostrará una diagonal perfecta, cuanto más se ciñan los puntos a esa diagonal perfecta más similar será distribución. El diagrama Q-Q consiste en comparar los cuantiles de la distribución observada con los cuantiles teóricos de una distribución normal con la misma media y desviación estándar que los datos. Cuanto más se aproximen los datos a una normal, más alineados están los puntos entorno a la recta.

Para la representación gráfica se puede hacer también uso de esta librería:

library(PASWR)
## Loading required package: e1071
##
## Attaching package: 'e1071'
## The following object is masked from 'package:Hmisc':
##
##     impute
EDA(P1)
## [1] "P1"

## Size (n)  Missing  Minimum   1st Qu     Mean   Median   TrMean   3rd Qu
##   19.000    0.000   75.000   92.000  105.579  111.000  105.579  118.000
##     Max.   Stdev.     Var.  SE Mean   I.Q.R.    Range Kurtosis Skewness
##  128.000   15.334  235.146    3.518   26.000   53.000   -1.064   -0.585
## SW p-val
##    0.042

Test de la varianza F de Fisher El F-test, también conocido como contraste de la razón de varianzas, contrasta la hipótesis nula de que dos poblaciones normales tienen la misma varianza. Es muy potente, detecta diferencias muy sutiles, pero es muy sensible a violaciones de la normalidad de las poblaciones

Se utiliza para saber si las varianzas de dos distribuciones son iguales. Como siempre: H0 las varianzas son iguales H1 las varianzas varían significativamente

var.test(P1,P2)
##
##  F test to compare two variances
##
## data:  P1 and P2
## F = 1.1874, num df = 18, denom df = 18, p-value = 0.7195
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.4574801 3.0820885
## sample estimates:
## ratio of variances
##           1.187432

Test de Bartlet Permite contrastar la igualdad de varianza en 2 o más poblaciones sin necesidad de que el tamaño de los grupos sea el mismo. Es más sensible que el test de Levene a la falta de normalidad, pero si se está seguro de que los datos provienen de una distribución normal, es la mejor opción.

data("iris")
a <- iris[iris$Species == "versicolor", "Petal.Length"]
b <- iris[iris$Species == "virginica", "Petal.Length"]
bartlett.test(list(a,b))
##
##  Bartlett test of homogeneity of variances
##
## data:  list(a, b)
## Bartlett's K-squared = 1.249, df = 1, p-value = 0.2637

Test de varianza Levene El test de Levene se caracteriza, además de por poder comparar 2 o más poblaciones, por permitir elegir entre diferentes estadísticos de centralidad :mediana (por defecto), media, media truncada. Esto es importante a la hora de contrastar la homocedasticidad dependiendo de si los grupos se distribuyen de forma normal o no.

library(car)
## Loading required package: carData
##
## Attaching package: 'carData'
## The following object is masked from 'package:PASWR':
##
##     Wool
data("iris")
iris <- dplyr::filter(.data = iris, Species %in% c("versicolor", "virginica"))
leveneTest(y = iris$Petal.Length, group = iris$Species, center = "median")
## Levene's Test for Homogeneity of Variance (center = "median")
##       Df F value Pr(>F)
## group  1  1.0674 0.3041
##       98

Este test se debe aplicar cuando se quieran comparar más de dos grupos como en un estudio ANOVA. Un p-value por debajo de 0.05 nos dice que las varianzas de los datos nos iguales, por tanto no se pueden usar test paramétricos como ANOVA.

library(car)
size<-c(25,22,28,24,26,24,22,21,23,25,26,30,25,24,21,27,28,23,25,24,20,22,24,23,22,24,20,19,21,22)
location<-c(rep("ForestA",10), rep("ForestB",10), rep("ForestC",10))
my.dataframe<-data.frame(size,location)
car::leveneTest(size~location, data=my.dataframe)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.4919 0.6168
##       27

Este test nos da un p-value muy alto lo que significa que no hay diferencias estadísticamente significativas. Se puede el parámetro center = mean para que use la media como centro en lugar de la mediana por defecto, que hace el test más robusto y menos sensible a outliers. En este caso el test se llama Brown-Forsythe

Test de varianza Fligner-Killeen

Si los datos no superasen el test de normalidad, deberíamos aplicar este otro método ya que es una prueba menos sensible a que los datos sean normales o cuando tengamos outliers no resueltos

#Debemos hacer también el test de Fligner ya que es una prueba menos sensible a que los datos sean normales
fligner.test(x = list(P1,P2))
##
##  Fligner-Killeen test of homogeneity of variances
##
## data:  list(P1, P2)
## Fligner-Killeen:med chi-squared = 0.04751, df = 1, p-value =
## 0.8275
#H0= Los datos tienen varianzas iguales H1: Los datos presentan varianzas distintas.

Datos emparejados

Lo último que podemos necesitar saber es si los sujetos de los grupos están o no emparejados. Si esto no lo sabemos a priori podemos hacer un análisis de correlación, si los datos tinen un alto nivel de correlación es un indicio de que estén emparejados.

cor(P1, P2)
## [1] 0.8838955
#H0= Los datos tienen varianzas iguales H1: Los datos presentan varianzas distintas.

Podemos comprobar gráficamente que están muy relacionadas

plot(P1, P2, pch=20, main = "P1 y P2")

a = 24 + 10*rbeta(150, 1.1, 1.1)  # generate fake data
b = 24 + 10*rbeta(150, 1.1, 1.1)
err = rnorm(150, 0, .1);  aa = a + err
plot(a, b, pch=20, main="Independent A and B");  plot(a, aa, pch=20, main="Paired A and AA")

Commo hemos podido comprobar, a y b están relaccionados a pesar de tener la misma distribución, mientras que P1 y P2 no parecen tener igual distribución aunque sí hemos visto que estaban relaccionadas:

qqplot(a,b,main="Q-QPlot A y B");qqplot(P1,P2,main="Q-QPlot P1 y P1")

Con toda esta infromación que hemos recabado, ya podemos confirmar si los datos son o no normales, si tienen o no distinta varianza y si se deben o no emparejar los resultados de los grupos o no.

U-Test o Wilcoxon Rank-Sum Test Datos NO EMPAREJADOS NO NORMALES Y VARIANZAS IGUALES

El test de Mann-Whitney-Wilcoxon (WMW), también conocido como Wilcoxon rank-sum test o u-test, es un test no paramétrico que contrasta si dos muestras proceden de poblaciones equidistribuidas.

La idea en la que se fundamenta este test es la siguiente: si las dos muestras comparadas proceden de la misma población, al juntar todas las observaciones y ordenarlas de menor a mayor, cabría esperar que las observaciones de una y otra muestra estuviesen intercaladas aleatoriamente. Por lo contrario, si una de las muestras pertenece a una población con valores mayores o menores que la otra población, al ordenar las observaciones, estas tenderán a agruparse de modo que las de una muestra queden por encima de las de la otra.

Al igual que ocurre con muchos test no paramétricos, el test de Mann-Whitney-Wilcoxon es menos potente que el t-test (tienen menos probabilidad de rechazar la H0 cuando realmente es falsa) ya que ignora valores extremos. En el caso de los t-test, al trabajar con medias, si los tienen en cuenta. Esto hace a su vez que el test de Mann-Whitney-Wilcoxon sea una prueba más robusta que los t-test. En concreto, la perdida de potencia es del 5%.

muestraX <- c( 1.1, 3.4, 4.3, 2.1, 7.0 , 2.5 )
muestraY <- c( 7.0, 8.0, 3.0, 5.0, 6.2 , 4.4 )
fligner.test(x = list(muestraX,muestraY))
##
##  Fligner-Killeen test of homogeneity of variances
##
## data:  list(muestraX, muestraY)
## Fligner-Killeen:med chi-squared = 0.07201, df = 1, p-value =
## 0.7884
#No hay evidencias en contra de la igualdad de varianzas.
wilcox.test(x = muestraX, y = muestraY, alternative = "two.sided", mu = 0,
            paired = FALSE, conf.int = 0.95)
## Warning in wilcox.test.default(x = muestraX, y = muestraY, alternative =
## "two.sided", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = muestraX, y = muestraY, alternative =
## "two.sided", : cannot compute exact confidence intervals with ties
##
##  Wilcoxon rank sum test with continuity correction
##
## data:  muestraX and muestraY
## W = 6.5, p-value = 0.07765
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  -4.9000274  0.4000139
## sample estimates:
## difference in location
##              -2.390035

En la salida devuelta por la función, el valor W equivale a U, con lo quee se puede obtener la probabilidad de superioridad como effect size calcullándolo de la siguiente manera:

p = U/N1*N2

en nuestro ejemplo p = 6,5/(6*6) = 0,18

La interpretación es que hay una probabilidad de un 18% de que un valor aleatorio del primer grupo sea mayor que un valor aleatorio del segundo grupo

Cuando hay ligaduras o ties(que al computar diferencias se repita mucho el mismo valor de la diferencia y no se puede rankear correctamente, por ejemplo de 30 resultados ordenados, al calcular la diferencia entre ellos hay muchos que se repite 2 de diferencia ), la función wilcox.test() no es capaz de calcular el p-value exacto por lo que devuelve una aproximación asumiendo que U se distribuye de forma ~ normal. En estos casos, o cuando los tamaños muestrales son mayores de 20 y se quiere la aproximación por la normal, es más recomendable emplear la función wilcox_test() del paquete coin.

require(coin)
## Loading required package: coin
# La función wilcox.test() del paquete coin requiere pasarle los argumentos en
# forma de función (~), por lo que los datos tienen que estar almacenados en
# forma de data frame.
datos <- data.frame(grupo = rep(c("A", "B"), c(6, 6)),
                    valores = c(muestraX, muestraY))
wilcox_test(valores ~ grupo, data = datos, distribution = "exact", conf.int=0.95)
##
##  Exact Wilcoxon-Mann-Whitney Test
##
## data:  valores by grupo (A, B)
## Z = -1.8447, p-value = 0.06926
## alternative hypothesis: true mu is not equal to 0
## 95 percent confidence interval:
##  -4.9  0.4
## sample estimates:
## difference in location
##                   -2.4
n1 <- length(muestraX)
n2 <- length(muestraY)
# La función devuelve el valor de Z, por lo que se puede calcular fácilmente el
# tamaño del efecto
tamanyo_efecto <- 1.845/sqrt(n1 + n2)
tamanyo_efecto
## [1] 0.5326056

Wilcoxon signed-rank test o Mann-Whitney-Wilcoxon (WMW) DATOS EMPAREJADOS NO NORMALES

El test no paramétrico prueba de los rangos con signo de Wilcoxon, también conocido como Wilcoxon signed-rank test, permite comparar poblaciones cuando sus distribuciones (normalmente interpretadas a partir de las muestras) no satisfacen las condiciones necesarias para otros test paramétricos. Es una alternativa al t-test de muestras dependientes cuando las muestras no siguen una distribución normal (muestran asimetría o colas) o cuando tienen un tamaño demasiado reducido para poder determinar si realmente proceden de poblaciones normales.

A la hora de elegir entre t-test o Wilcoxon signed-rank test, es importante tener en cuenta que, el problema de las muestras pequeñas, no se soluciona con ninguno de los dos. Si el tamaño de las muestras es pequeño, también lo es la calidad de la inferencia que se puede hacer con ellas. Ahora bien, existen dos situaciones en las que, a priori, se puede recomendar utilizar un Wilcoxon signed-rank test antes que un t-test:

Si el tamaño de las muestras es suficientemente grande para determinar (por métodos gráficos o contrastes de hipótesis) que la distribución de las poblaciones a comparar no es de tipo normal, en tal caso, los t-test no son adecuados, por lo que mejor emplear un Wilcoxon signed-rank test (Bootstrapping, regresión cuantílica, o test de permutación también podrían ser otras alternativa).

Si el tamaño de las muestras no permite determinar con seguridad si las poblaciones de las que proceden se distribuyen de forma normal, y no se dispone de información que pueda orientar sobre la naturaleza de las poblaciones de origen (estudios anteriores, que sea un tipo de variable que se sabe que se distribuye casi siempre de forma normal???), entonces es más apropiado el Wilcoxon signed-rank test ya que no requiere asumir la normalidad de las poblaciones.

El test Wilcoxon signed-rank test se caracteriza por:

Es frecuente encontrar descrito que, el Wilcoxon signed-rank test, compara la mediana de las diferencias, sin embargo, esto solo es correcto bajo determinadas condiciones. A modo general, el Wilcoxon signed-rank test compara si las diferencias entre pares de datos siguen una distribución simétrica entorno a un valor. Si dos muestras proceden de la misma población, es de esperar que las diferencias entre cada par de observaciones se distribuyan de forma simétrica entorno al cero.

Trabajan sobre rangos de orden, es decir, utilizan las posiciones que ocupan los datos una vez ordenados. Por lo tanto, solo es aplicable a variables cuyos valores se pueden ordenar.

Tienen menos poder estadístico (menor probabilidad de rechazar la hipótesis nula cuando realmente es falsa) ya que ignoran valores extremos. En el caso de los t-test, al trabajar con medias, sí los tienen en cuenta. Esto a su vez, hace que el Wilcoxon signed-rank test sean una prueba más robusta que el t-test.

antes       <- c( 2, 5, 4, 6, 1, 3 )
despues     <- c( 5, 6, 2, 7, 1, 6 )

wilcox.test(x = antes, y = despues, alternative = "two.sided", mu = 0, paired = TRUE)
## Warning in wilcox.test.default(x = antes, y = despues, alternative =
## "two.sided", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = antes, y = despues, alternative =
## "two.sided", : cannot compute exact p-value with zeroes
##
##  Wilcoxon signed rank test with continuity correction
##
## data:  antes and despues
## V = 3, p-value = 0.2763
## alternative hypothesis: true location shift is not equal to 0

En la salida devuelta por wilcox.test(), el estadístico se denomina V en lugar de W.

p = 3/(6*6) = 0,08

Cuando hay ligaduras o ties, wilcox.test() no es capaz de calcular el p-value exacto, por lo que devuelve un p-value aproximado asumiendo que W se distribuye de forma aproximadamente normal. En estos casos, o cuando los tamaños muestrales son mayores de 25, es recomendable emplear la función wilcoxsigned_test() del paquete coin, que devuelve el valor exacto de p-value en lugar de una aproximación.

require(coin)
# La función wilcoxsigned_test() del paquete coin requiere pasarle los
# argumentos en forma de función (~), por lo que los datos tienen que estar
# almacenados en forma de data frame.

datos <- data.frame(antes = antes, despues = despues)
wilcoxsign_test(antes ~ despues, data = datos, distribution = "exact")
##
##  Exact Wilcoxon-Pratt Signed-Rank Test
##
## data:  y by x (pos, neg)
##   stratified by block
## Z = -1.272, p-value = 0.25
## alternative hypothesis: true mu is not equal to 0

WILCOXON VS T-TEST

Suele recomendarse emplear el test de Mann???Whitney???Wilcoxon en lugar del t-test cuando los tamaños muestrales son pequeños y no se tiene evidencias de que las poblaciones de origen siguen una distribución normal. Si bien esta práctica está bastante fundamentada, no hay que confundirla con la de utilizar el test de Mann???Whitney???Wilcoxon como alternativa al t-test siempre que no se cumpla la normalidad y sin tener en cuenta el tamaño muestral. A medida que el número de observaciones aumenta, también lo hace la robustez del t-test frente a desviaciones de la normalidad. Por otro lado, el test de Mann???Whitney???Wilcoxon incrementa su sensibilidad a diferencias más allá de las medianas. Por ejemplo, aumenta el poder estadístico de detectar diferencias significativas en la probabilidad de que observaciones de un grupo superen a las del otro debido únicamente a diferencias en la dispersión de las poblaciones de origen y no porque sus medianas sean distintas. Si el objetivo del estudio es identificar cualquier diferencia distribucional, esto no supone un problema, pero si lo que se quiere comparar son las medianas, se obtendrán p-values que no se corresponden con la pregunta que el investigador quiere responder.

Considérese dos muestras que proceden de dos poblaciones distintas, una con distribución normal y otra con distribución log-normal (asimetría derecha). Ambas poblaciones tienen la misma mediana pero claramente siguen distribuciones distintas.

x <- seq(0, 5, length = 1000)
y <- dlnorm(x = x, meanlog =  0.5, sdlog = 1.7)
plot(x, y, type = "l", lty = 1, xlab = "x", col = "blue", ylab = "Densidad",
     main = "Distribuciones lognormal y normal con misma mediana", xlim=c(-5, 5))
x_2 <- seq(-5, 5, length = 1000)
y_2 <- dnorm(x = x_2, mean = 1.648721, sd = 1)
lines(x_2, y_2, col = "red")
legend("topright",
       legend = c("mediana lognormal = 0.5", "mediana normal = 0.5"),
       col = c("blue", "red"), lty = 1,cex = 0.7)

Supóngase ahora que el investigador, que desconoce la distribución real de las poblaciones, obtiene dos muestras de 20 observaciones cada una.

set.seed(888)
muestra_a <- rlnorm(n = 20, meanlog =  0.5, sdlog = 1.7)
muestra_b <- rnorm(n = 20, mean = 1.648721, sd = 1)
par(mfrow = c(1,2))
hist(muestra_a, breaks = 10, main = "", col = "blue")
hist(muestra_b, breaks = 10, main = "", col = "firebrick")

En base a la falta de normalidad de la muestra a, el investigador considera que el t-test no es adecuado para comparar la localización de las poblaciones y decide que quiere comparar medianas con el test de Mann-Whitney-Wilcoxon.

wilcox.test(muestra_a, muestra_b, paired = FALSE)
##
##  Wilcoxon rank sum test
##
## data:  muestra_a and muestra_b
## W = 196, p-value = 0.9254
## alternative hypothesis: true location shift is not equal to 0

El resultado le indica que no hay evidencias para considerar que la localización de las poblaciones es distinta y concluye que las medianas de ambas poblaciones son iguales.

Véase lo que ocurre si el tamaño muestral se incrementa de 20 a 1000 observaciones.

set.seed(888)
muestra_a <- rlnorm(n = 1000, meanlog = 0.5, sdlog = 1.7)
muestra_b <- rnorm(n = 1000, mean = 1.648721, sd = 1)
wilcox.test(muestra_a, muestra_b, paired = FALSE)
##
##  Wilcoxon rank sum test with continuity correction
##
## data:  muestra_a and muestra_b
## W = 551900, p-value = 5.842e-05
## alternative hypothesis: true location shift is not equal to 0

En este caso, la conclusión obtenida es totalmente opuesta, hay muchas evidencias en contra de la hipótesis nula de que ambas poblaciones tienen la misma localización por lo que el investigador podría concluir erróneamente que las medianas de ambas poblaciones son distintas.

Este ejemplo pone de manifiesto la importancia que tiene el tamaño muestral en la potencia del test de Mann???Whitney???Wilcoxon y la necesidad de entender que, si las distribuciones no son iguales a excepción de su localización, el test es muy sensible detectando diferencias en las distribuciones a pesar de que sus medianas sean las mismas. En este caso, el investigador solo podría concluir que hay evidencias de que la distribución de las dos poblaciones se diferencia en algún aspecto, pero no podría concretar cuál.

CONTRASTE DE HIPÓTESIS DATOS CATEGÓRICOS

Los contrastes de hipótesis para variables cualitativas se realizan mediante test de frecuencia o proporciones. Dentro de esta categoría existen distintos tipos de test, la utilización de uno u otro depende de qué tipo de información se quiera obtener:

Test de distribución esperada o goodness of fit: Se emplean para comparar la distribución observada frente a una distribución esperada o teórica. Test de diferencia de frecuencias o test de independencia: Se emplean para estudiar si la frecuencia de observaciones es significativamente distinta entre dos o más grupos. En los test de goodness of fit solo hay una variable asociada a cada observación, mientras que en los test de independencia hay dos variables asociadas a cada observación. También se emplean distintos test dependiendo del tipo de datos (independientes o dependientes) con los que se vaya a trabajar. Las siguientes tablas muestran algunos de los más empleados.

Tabla de métodos Los test exactos calculan la probabilidad de obtener los resultados observados de forma directa generando todos los posibles escenarios y calculando la proporción en los que se cumple la condición estudiada (son test de permutaciones). Los test aproximados calculan primero un estadístico y luego emplean la distribución teórica de dicho estadístico para obtener la probabilidad de que adquiera valores iguales o más extremos.

Existe bastante controversia en cuanto a si se deben de utilizar test exactos o aproximados. En la era pre-computacional, los test exactos se complicaban mucho cuando el tamaño total de muestras aumentaba, sin embargo, por medio de la computación esta barrera se ha eliminado. Los test exactos son más precisos cuando el tamaño total de observaciones es bajo o alguno de los grupos tiene pocas observaciones, una vez alcanzado un número alto de observaciones las diferencias son mínimas. En el libro Handbook of Biological Statistics John H. McDonald se recomienda utilizar test exactos cuando el número total de observaciones es menor a 1000 o cuando, aunque el número total sea mayor a 1000, haya algún grupo cuyo número de eventos esperados sea pequeño (normalmente menor que 5). En el caso de aplicar test aproximados sobre tamaños pequeños se suelen emplear correcciones, las más frecuentes son la corrección de continuidad de Yate o la corrección de William.

Se puede considerar a los test basados el estadístico ??2 como una generalización del contraste de proporciones basado en la aproximación a la normal (Z-test de una proporción y Z-test de dos proporciones) cuando hay 2 o más variables cualitativas o alguna de ellas tiene 2 o más niveles. En aquellos casos en los que ambos test se pueden aplicar, el resultado de un Z-test y un test ??2 es equivalente. Esto es debido a que en la distribución chi-cuadrado con 1 grado de libertad el estadístico ??2 es igual al estadístico Z de una distribución normal, elevado al cuadrado. ##1 DISTRIBUCIÓN ESPERADA## ###TEST BINOMIAL EXACTO### El test binomial exacto se emplea para estudiar si la proporción de eventos verdaderos de una variable de tipo binomial se diferencia significativamente de la frecuencia teórica con la que se esperaría que apareciesen. En R puede realizarse un test binomial exacto mediante la función binom.test(), obteniendo tanto el p-value como el intervalo de confianza del valor de p en la población.

Un estudio sugiere que en una ciudad están naciendo más hombres que mujeres. Para determinar si es cierto se selecciona una muestra aleatoria de los niños nacidos durante los últimos 2 años y se identifica para cada uno el género. ¿Existen diferencias significativas para un nivel de significancia del 5%?

datos <- c("masculino", "masculino", "masculino", "masculino", "femenino",
           "masculino", "masculino", "masculino", "femenino", "femenino",
           "femenino", "femenino", "femenino", "femenino", "masculino",
           "masculino", "femenino", "masculino", "femenino", "masculino",
           "femenino", "masculino", "masculino", "masculino", "femenino",
           "masculino", "masculino", "femenino", "masculino", "femenino")
tabla <- table(datos)
tabla
## datos
##  femenino masculino
##        13        17
binom.test(x = tabla, alternative = "two.sided", conf.level = 0.95)
##
##  Exact binomial test
##
## data:  tabla
## number of successes = 13, number of trials = 30, p-value = 0.5847
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.2546075 0.6257265
## sample estimates:
## probability of success
##              0.4333333

El test binomial no encuentra diferencias significativas entre las observaciones y lo que cabría esperar si la verdadera probabilidad de nacimiento es del 50% para ambos sexos.

Más de 2 niveles

En los estudios genéticos de Mendel se establece que la descendencia de un determinado cruce de plantas da lugar a cuatro tipos de semillas (amarillas-lisas, verdes-lisas, amarillas-rugosas y verdes-rugosas) con una proporción teórica de 9:3:3:1. Para demostrar si es cierto se recoge una muestra de 556 semillas y se observa que las cantidades para cada grupo son de 315, 108, 101, 32 respectivamente. ¿Se corresponden las observaciones con el modelo esperado?

observaciones <- c(315, 108, 101, 32)
freq_teoricas <- c(9,3,3,1)

require(XNomial)
## Loading required package: XNomial
XNomial::xmulti(obs = observaciones, expr = freq_teoricas, statName = "Prob") 
##
## P value (Prob) = 0.9382

La función xmulti() también calcula el p-value a partir de modelos aproximados (explicados más adelante) y crea una representación gráfica que permite ver como de bien se ajusta la distribución chi-cuadrado a las probabilidades exactas.

xmulti(obs = observaciones, expr = freq_teoricas, detail = 3, histobins = TRUE)
##
## P value  (LLR)  =  0.9261
## P value (Prob)  =  0.9382
## P value (Chisq) =  0.9272
##
## Observed:  315 108 101 32
## Expected ratio:  9 3 3 1
## Total number of tables:  28956759

El p-value obtenido indica que las observaciones se asemejan mucho a las esperadas según la distribución de frecuencias del modelo teórico (Ho). No hay evidencias para rechazar la hipótesis nula.

CHI2 GOODNESS OF FIT

El test chi2 goodness of fit es la alternativa aproximada al test binomial (dos niveles) o multinomial (>2 niveles) normalmente utilizada cuando el número de observaciones o de grupos es demasiado alto para poder emplear los test exactos. Permite comparar la distribución de las observaciones en los diferentes niveles de una variable cualitativa con lo esperado acorde a una distribución hipotética Ho. Se trata por lo tanto de una expansión del Z-test para una proporción cuando la variable estudiada tiene dos o más niveles. Cuando la variable cualitativa tiene dos niveles, el test chi2 goodness of fit y el Z-test para una proporción son equivalentes.

El funcionamiento del test consiste en cuantificar y sumar las diferencias entre el número de eventos observadas en cada nivel con respecto al número esperado acorde a H0. El valor de este sumatorio se asigna a un estadístico llamado chi2 y se emplea su relación matemática con la distribución chi2 para estimar la probabilidad de obtener un valor igual o más extremo.

H0: Las observaciones siguen una distribución acorde a la distribución del modelo de referencia. Las desviaciones son debidas al azar.

Ha: Las observaciones no siguen una distribución acorde a la distribución del modelo de referencia. Las desviaciones son mayores a las esperadas por azar.

Un estudio intenta comparar si existe discriminación racial en la formación de jurados populares. Para ello se compara la distribución racial en el país con la distribución de los miembros de los últimos 500 jurados.

Ejemplo

Ejemplo

chisq.test(x = c(1920, 347, 19, 84, 130), p = c(80.29, 12.06, 0.79, 2.92, 3.94) / 100)
##
##  Chi-squared test for given probabilities
##
## data:  c(1920, 347, 19, 84, 130)
## X-squared = 22.419, df = 4, p-value = 0.0001654

En caso de que no se cumpla el número mínimo de eventos por nivel, se recurre a la simulación:

chisq.test(x = c(1920, 347, 19, 84, 130),
           p = c(80.29, 12.06, 0.79, 2.92, 3.94)/100,
           simulate.p.value = TRUE, B = 5000)
##
##  Chi-squared test for given probabilities with simulated p-value
##  (based on 5000 replicates)
##
## data:  c(1920, 347, 19, 84, 130)
## X-squared = 22.419, df = NA, p-value = 0.0005999

Además del valor del p-value se recomienda incluir una representación gráfica de los valores observados y esperados.

require(ggplot2)
datos <- data.frame(raza=rep(c("blanca","negra","latina","asiática","otra"), 2),
                    tipo_valor = rep(c("observado", "esperado"), each = 5),
                    eventos = c(1920, 347, 19, 84, 130, 2007.25, 301.50,
                                19.75, 73.00, 98.50))
ggplot(data = datos, aes(x = raza, y = eventos)) +
geom_bar(stat = "identity", aes(fill = tipo_valor), position = "dodge") +
theme_bw()

Se confirma tanto por un pvalue muy bajo como gráficamente que no se está siguiendo la proporción esperada

2 COMPARACIÓN DE GRUPOS

TEST EXACTO DE FISHER

La prueba de Fisher es el test exacto utilizado cuando se quiere estudiar si existe asociación entre dos variables cualitativas, es decir, si las proporciones de una variable son diferentes dependiendo del valor que adquiera la otra variable. En la gran mayoría de casos, el test de Fisher se aplica para comparar dos variables categóricas con dos niveles cada una (tabla 2x2). Es posible utilizarlo con tablas 2xK niveles pero los requerimientos de cálculo son altos.

Ho: Las variables son independientes por lo que una variable no varía entre los distintos niveles de la otra variable.

Ha: Las variables son dependientes, una variable varía entre los distintos niveles de la otra variable.

El test de Fisher es más preciso que sus equivalentes aproximados (test chi-square de independencia o G???test de independencia) cuando el número de eventos esperado por nivel es pequeño. Se recomienda utilizarlo siempre que sea posible (tiempo de computación) aunque para observaciones totales >1000 los resultados de los test aproximados son muy parecidos.

Es importante tener en cuenta que el test de Fisher está diseñado para situaciones en las que las frecuencias marginales de filas y columnas (los totales de cada fila y columna) son fijas, se conocen de antemano. Esta condición es relevante en los experimentos biológicos ya que no es común poder cumplirla. Si esta condición no se satisface el test de Fisher deja de ser exacto, por lo general pasando a ser más conservativo. En varios artículos se menciona que el test de Barnard es más potente que el de Fisher cuando las frecuencias marginales no son fijas. También parece ser que aunque el test deja de ser exacto no significa que no se pueda aplicar.

Ejemplo de experimentos con y sin frecuencias marginales fijas:

Frecuencias marginales fijas: Supóngase que se quiere saber si la preferencia que tienen dos especies de pájaros (estorninos y gorriones) para refugiarse en casetas artificiales es diferente dependiendo del material de fabricación (madera o metal). Para ellos se disponen en una pajarera 5 casetas de metal y 5 de madera y se sueltan en el interior de la jaula 4 gorriones y 6 estorninos. En este experimento se sabe que las frecuencias marginales van a ser 5, 5, 4, 6 lo que no se sabe es como se van a distribuir las observaciones dentro de la tabla.

Tabla de métodos

Tabla de métodos

Dado que el test de Fisher contrasta si las variables están relacionadas, al tamaño del efecto se le conoce como fuerza de asociación. Existen múltiples medidas de asociación, entre las que destacan phi o Cramer???s V. Los límites empleados para su clasificación son:

pequeño: 0.1 mediano: 0.3 grande: 0.5 En R se pueden calcular mediante la función assocstats() del paquete vcd.

Se quiere estudiar si la reacción alérgica a un compuesto y una determinada mutación en un gen están relacionados. Para ello se realiza un test alérgico sobre un grupo de individuos seleccionados al azar y se genotipa el estado del gen de interés ¿Existe un diferencia significativa en la incidencia de la mutación entre los alérgicos y no alérgicos?

datos <- data.frame( sujeto = c("No alérgico", "No alérgico", "No alérgico",
                                "No alérgico", "alérgico", "No alérgico",
                                "No alérgico", "alérgico", "alérgico",
                                "No alérgico", "alérgico", "alérgico",
                                "alérgico", "alérgico", "alérgico",
                                "No alérgico", "No alérgico", "No alérgico",
                                "No alérgico", "alérgico", "alérgico",
                                "alérgico", "alérgico", "No alérgico",
                                "alérgico", "No alérgico", "No alérgico",
                                "alérgico","alérgico", "alérgico"),
                     mutacion = c(FALSE, FALSE, FALSE, FALSE, TRUE, FALSE,
                                  FALSE, FALSE, TRUE, TRUE, TRUE, TRUE, TRUE,
                                  TRUE, FALSE,  FALSE, TRUE, FALSE, TRUE, FALSE,
                                  TRUE, FALSE,FALSE, FALSE, TRUE, FALSE, FALSE,
                                  TRUE, FALSE, TRUE))

paste("El test de fisher trabaja con tablas")
## [1] "El test de fisher trabaja con tablas"
tabla <- table(datos$sujeto, datos$mutacion, dnn = c("Sujeto", "Estado gen"))
tabla
##              Estado gen
## Sujeto        FALSE TRUE
##   alérgico        6   10
##   No alérgico    11    3
fisher.test(x = tabla, alternative = "two.sided")
##
##  Fisher's Exact Test for Count Data
##
## data:  tabla
## p-value = 0.03293
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.02195148 1.03427479
## sample estimates:
## odds ratio
##  0.1749975
library(vcd)
## Loading required package: grid
paste0("Fuerza de asociación")
## [1] "Fuerza de asociación"
assocstats(x = tabla)
##                     X^2 df P(> X^2)
## Likelihood Ratio 5.3356  1 0.020894
## Pearson          5.1293  1 0.023525
##
## Phi-Coefficient   : 0.413
## Contingency Coeff.: 0.382
## Cramer's V        : 0.413

En este ejemplo no se satisface la condición de frecuencias marginales fijas y por lo tanto el test de Fisher no es exacto. Aun así, hay evidencias para rechazar la H0 y considerar que las dos variables sí están relacionadas. El tamaño de la fuerza de asociación (tamaño de efecto) cuantificado por phi o Cramers V es mediano.

Chi cuadrado de Pearson

El test chi2 de independencia, también conocido como chi2 de Pearson se emplea para estudiar si existe asociación entre dos variables categóricas, es decir, si las proporciones de una variable son diferentes dependiendo del valor que adquiera la otra variable, cuando los datos son independientes. Se trata por lo tanto de una expansión del Z-test para dos proporcines cuando una de las variables estudiadas tiene dos o más niveles. Cuando ambas variables tienen dos niveles (tabla 2x2) ambos test chi2 goodness of fit y Z-test para dos proporción son equivalentes.

Es el test aproximado equivalente a su versión exacta test de Fisher. Debido a los requerimientos de cálculo del test de Fisher, cuando hay muchas observaciones o muchos niveles, se emplea el test ??2 de independencia. Es importante tener en cuenta que cuando el número de observaciones esperadas para alguno de los niveles es igual o menor a 5 la aproximación por el test ??2 no es buena.

El test de independencia cuantifica y sumariza cómo de distinto es el número de eventos observados en cada nivel con respecto al número esperado acorde con Ho. Esto permite identificar si la desviación total es mayor que la que cabría esperar simplemente por azar.

H0: Las variables son independientes por lo que una variable no varía entre los distintos niveles de la otra variable.

Ha: Las variables son dependientes, una variable varía entre los distintos niveles de la otra variable.

ndependencia: las observaciones de la muestra deben ser independientes unas de otras.

Muestreo aleatorio. Tamaño de la muestra < 10% población. Cada observación contribuye únicamente a uno de los niveles. Tamaño: cada nivel debe tener al menos 5 eventos esperados (acorde a Ho) y el número de observaciones totales (n) >30. En caso de no cumplirse esta condición, el test pierde precisión.

En el libro Bioestadística de Francisxa Ruis Díaz consideran que esta regla es muy estricta y rara vez se cumple en la práctica. Proponen unas condiciones más relajadas con las que no se pierde demasiada precisión: para ningún nivel el número de eventos esperados acorde a H0 es menor de 1 y como máximo un 20% de los niveles tiene menos de 5 eventos esperados.

En caso de no cumplirse esta condición o estar en el límite, se recurre a los test exactos o, si no es posible, a la simulación. Aun cuando se cumplen las condiciones, son más precisos los test exactos y por lo tanto más recomendables.

fila1 <- c(81, 103, 147)
fila2 <- c(359, 326, 277)
tabla <- as.table(rbind(fila1, fila2))
dimnames(tabla) = list(Peso = c("Obeso","No obeso"),
                       Estado_civil = c("soltero","pareja","casado"))
tabla
##           Estado_civil
## Peso       soltero pareja casado
##   Obeso         81    103    147
##   No obeso     359    326    277
chisq.test(x = tabla)
##
##  Pearson's Chi-squared test
##
## data:  tabla
## X-squared = 30.829, df = 2, p-value = 2.021e-07
paste("Solución mediante simulación")
## [1] "Solución mediante simulación"
chisq.test(tabla, simulate.p.value = TRUE, B = 5000)
##
##  Pearson's Chi-squared test with simulated p-value (based on 5000
##  replicates)
##
## data:  tabla
## X-squared = 30.829, df = NA, p-value = 2e-04

Dado que el test ha resultado significativo, se quiere determinar para que niveles el número de observaciones difiere significativamente de lo esperado.

chisq.test(x = tabla)$residuals
##           Estado_civil
## Peso          soltero     pareja     casado
##   Obeso    -2.9809729 -0.6509186  3.6914422
##   No obeso  1.7485759  0.3818151 -2.1653222
chisq.test(x = tabla)$stdres
##           Estado_civil
## Peso          soltero     pareja     casado
##   Obeso    -4.2549504 -0.9231681  5.2203206
##   No obeso  4.2549504  0.9231681 -5.2203206

Las mayores desviaciones respecto a los valores esperados se dan en los niveles de soltero y casado. Para tener un estudio más exacto de si hay diferencias significativas en estos niveles, se divide la tabla en tablas 2X2 y se repite el test ??2 corrigiendo el nivel de significancia.

solteros_casados <- tabla[, c(1,3)]
solteros_casados
##           Estado_civil
## Peso       soltero casado
##   Obeso         81    147
##   No obeso     359    277
chisq.test(solteros_casados)
##
##  Pearson's Chi-squared test with Yates' continuity correction
##
## data:  solteros_casados
## X-squared = 28.56, df = 1, p-value = 9.083e-08
assocstats(solteros_casados)
##                     X^2 df   P(> X^2)
## Likelihood Ratio 29.687  1 5.0775e-08
## Pearson          29.391  1 5.9140e-08
##
## Phi-Coefficient   : 0.184
## Contingency Coeff.: 0.181
## Cramer's V        : 0.184

Se confirma que al menos entre los grupos casados y solteros sí hay una asociación significativa de las variables estado civil y obesidad con un tamaño de asociación Cramers V pequeño. Por lo tanto se puede afirmar que el porcentaje de gente obesa está asociado al estado civil.

Test de McNemar DATOS PAREADOS

El test de McNemar es la alternativa a los test ??2 de Pearson y al test exacto de Fisher cuando: los datos son pareados, se trata de tablas 2x2 y ambas variables son dicotómicas (binomiales). El test de McNemar estudia si la probabilidad de evento verdadero para una variable es igual en los dos niveles de otra variable.

Tabla de McNemar

Tabla de McNemar

Supóngase que se quiere comprobar si un tratamiento de hipnosis es capaz de hacer que las personas contesten ???Sí??? con mayor frecuencia. Para ello se selecciona un grupo de individuos a los que se les realiza una pregunta cuya respuesta puede ser SI/NO antes y después de someterse al tratamiento de hipnosis.

Ho: La respuesta de los sujetos es independiente del tratamiento. La proporción de sujetos que pasan de responder Sí a No es igual a la proporción de sujetos que pasan de responder No a Sí.

Ha: El tratamiento sí influye en la respuesta por lo que la proporción de sujetos que pasan de No a Sí es diferente de la de sujetos que pasan de Sí a No.

datos <- data.frame( sujeto = rep(1:15, each = 2),
                     tratamiento = c("pre" , "post" , "pre" , "post" , "pre" ,
                                     "post" , "pre" , "post" , "pre" , "post" ,
                                     "pre" , "post" , "pre" , "post" , "pre" ,
                                     "post" , "pre" , "post" , "pre" , "post" ,
                                     "pre" , "post" , "pre" , "post" , "pre" ,
                                     "post" , "pre" , "post" , "pre" , "post"),
                     respuesta = c("NO", "SI", "SI", "SI", "NO", "SI", "SI",
                                   "NO", "SI", "SI", "NO", "SI", "NO", "SI",
                                   "NO", "SI", "NO", "SI", "SI", "SI", "NO",
                                   "NO", "SI", "SI", "NO", "SI", "NO", "NO",
                                   "NO", "SI"))
library(tidyr)
datos <- spread(data = datos, key = tratamiento, value = respuesta)
head(datos)
##   sujeto post pre
## 1      1   SI  NO
## 2      2   SI  SI
## 3      3   SI  NO
## 4      4   NO  SI
## 5      5   SI  SI
## 6      6   SI  NO
tabla <- table(Pre_Tratamiento = datos$pre, Post_Tratamiento = datos$post)
tabla
##                Post_Tratamiento
## Pre_Tratamiento NO SI
##              NO  2  8
##              SI  1  4
mcnemar.test(tabla)
##
##  McNemar's Chi-squared test with continuity correction
##
## data:  tabla
## McNemar's chi-squared = 4, df = 1, p-value = 0.0455
binom.test(x = 1, n = 1 + 8, p = 0.5 )
##
##  Exact binomial test
##
## data:  1 and 1 + 8
## number of successes = 1, number of trials = 9, p-value = 0.03906
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.002809137 0.482496515
## sample estimates:
## probability of success
##              0.1111111

Existen evidencias significativas para rechazar Ho en favor de que sí existe relación entre el tratamiento y la respuesta de los sujetos.

Test de Cochran DATOS PAREADOS MÁS DE 2 NIVELES

El test Q-Cochran es el equivalente al test de McNemar para más de dos grupos. Permite estudiar la independencia entre varias muestras pareadas, es decir, si la distribución de una variable binomial es la misma en todos los grupos. Un test Q-Conchran de dos grupos es equivalente a un test de McNemar.

H0: La proporción de eventos verdaderos es la misma en los diferentes grupos. Ha: La proporción de eventos verdaderos es distinta en los diferentes grupos.

Si el test Q-Cochran es significativo, implica que al menos dos grupos difieren entre ellos. Para identificar cuales son se recurre a comparar dos a dos los diferentes grupos mediante el test de McNemar haciendo corrección de significancia (Bonferroni, Holm u otra). Esto se llama análisis Post-Hoc

No existe un método general de calcular la fuerza de asociación para el Q-Cochran test. Lo que se hace en su lugar es calcular la fuerza de asociación para cada comparación en el análisis post hoc, en este caso el de los test McNemar.

Respuesta <- c(1,0,1,0,0,1,0,1,0,0,1,1,0,0,1,1,1,1,0,0,1,0,1,1,0,1,1,0,1,1)
Sujeto <- factor(c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6,7,7,7,8,8,8,9,9,9,10,
                   10,10))
Canal <- factor(rep(c("inicial","anuncio","internet"), 10))
datos <- data.frame(Sujeto,Canal,Respuesta)

Se trata de una variable binomial cuya distribución (proporción de 1???s y 0???s) se quiere estudiar en más de dos grupos. Dado que los datos están pareados se emplea un test Q-Cochran.

 require(coin)
symmetry_test(Respuesta ~ factor(Canal) | factor(Sujeto), data = datos,teststat = "quad") 
##
##  Asymptotic General Symmetry Test
##
## data:  Respuesta by
##   factor(Canal) (anuncio, inicial, internet)
##   stratified by factor(Sujeto)
## chi-squared = 8.2222, df = 2, p-value = 0.01639

Existen evidencias significativas de que la intención de compra es distinta al menos entre dos grupos.

Comparaciones Post-Hoc

Se comparan dos a dos las 3 categorías del modo de información mediante test de McNemar.

 require(coin)
symmetry_test(Respuesta ~ factor(Canal) | factor(Sujeto), data = datos,teststat = "quad") 
##
##  Asymptotic General Symmetry Test
##
## data:  Respuesta by
##   factor(Canal) (anuncio, inicial, internet)
##   stratified by factor(Sujeto)
## chi-squared = 8.2222, df = 2, p-value = 0.01639
require(coin)
require(dplyr)
## Loading required package: dplyr
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
##
##     recode
## The following objects are masked from 'package:Hmisc':
##
##     src, summarize
## The following object is masked from 'package:MASS':
##
##     select
## The following objects are masked from 'package:stats':
##
##     filter, lag
## The following objects are masked from 'package:base':
##
##     intersect, setdiff, setequal, union
require(tidyr)
# Inicial vs anuncio
inicial_anuncio <- datos %>% filter(Canal %in% c("inicial","anuncio")) %>% 
                   spread(key = Canal, value = Respuesta)
tabla_inicial_anuncio <- table(inicial = inicial_anuncio$inicial,
                               anuncio = inicial_anuncio$anuncio)
tabla_inicial_anuncio
##        anuncio
## inicial 0 1
##       0 3 5
##       1 1 1
p_value_12 <- mcnemar.test(tabla_inicial_anuncio)
p_value_12
##
##  McNemar's Chi-squared test with continuity correction
##
## data:  tabla_inicial_anuncio
## McNemar's chi-squared = 1.5, df = 1, p-value = 0.2207
# Inicial vs internet
inicial_internet <- datos %>% filter(Canal %in% c("inicial","internet")) %>% 
                    spread( key = Canal, value = Respuesta)
tabla_inicial_internet <- table(inicial = inicial_internet$inicial,
                                internet = inicial_internet$internet)
tabla_inicial_internet
##        internet
## inicial 0 1
##       0 1 7
##       1 0 2
p_value_13 <- mcnemar.test(tabla_inicial_internet)
p_value_13
##
##  McNemar's Chi-squared test with continuity correction
##
## data:  tabla_inicial_internet
## McNemar's chi-squared = 5.1429, df = 1, p-value = 0.02334
# anuncio vs internet
anuncio_internet <- datos %>% filter(Canal %in% c("anuncio","internet")) %>%
                    spread(key = Canal, value = Respuesta)
tabla_anuncio_internet <- table(anuncio = anuncio_internet$anuncio,
                                internet = anuncio_internet$internet)
tabla_anuncio_internet
##        internet
## anuncio 0 1
##       0 0 4
##       1 1 5
p_value_23 <- mcnemar.test(tabla_anuncio_internet)
p_value_23
##
##  McNemar's Chi-squared test with continuity correction
##
## data:  tabla_anuncio_internet
## McNemar's chi-squared = 0.8, df = 1, p-value = 0.3711

Se corrige el nivel de significancia, en este caso se selecciona el método de holm.

p.adjust(c(p_value_12$p.value, p_value_13$p.value, p_value_23$p.value), method = "holm")
## [1] 0.44134272 0.07002661 0.44134272

A pesar de que el test Q-Cochran ha resultado significativo, las comparaciones dos a dos con corrección de holm no encuentran ninguna diferencia significativa. Es importante mencionar que en este caso no se satisfacían las condiciones para un test de McNemar por lo que habría que utilizar para las comparaciones post-hoc el test binomial.

P-VALUE

Hemos visto cómo calcular comparar las media de dos grupos y decir cuál es su diferencia en un 95% de los casos con un inervalo de confianza, y hemos visto como contrastar los resultados de dos grupos, comporbando si el valor de uno estaba o no dentro del 5% de posibilidades del otro para poder descartar una hipótesis frente a la otra. Estos dos conceptos los unimos para porporcionar en un contraste de hipótesis el p-value, con el que compararmos nuestra hipótesis con la posibilidad de que al azar obtengamos ese mismo valor, proporcionando un intervalo de confianza.

Nos sirve para dictaminar que el resultado obtenido es estadísticamnete significativo.

qownnotes-media-AeDUXI

qownnotes-media-AeDUXI

Es la probabilidad de obtener un resultado al menos tan extremo como el que realmente se ha obtenido (valor del estadístico calculado), suponiendo que la hipótesis nula es cierta. EL pvalue por tanto es una medida de significación estadística (un resultado o efecto es estadísticamente significativo cuando es improbable que haya sido debido al azar).

Una característica de p es que su valor depende del tamaño de la muestra a partir de la cual se determina: a mayor tamaño muestral mayor es la probabilidad de encontrar diferencias significativas entre los grupos que se comparan. De ahí la importancia de que en todo estudio deba realizarse el cálculo del tamaño muestral necesario en base a las caractrísticas de dicho estudio con anterioridad a la recogida y el análisis de los datos.

qownnotes-media-WX4036

qownnotes-media-WX4036

Todos los libros de estadística tienen una tabla detrás con los pvalues asociado a cada uno de los valores de Zscore de distintas distribuciones, en este caso distribución normal

qownnotes-media-Tt4036

qownnotes-media-Tt4036

En primer paso por tanto para probar nuestra hipótesis será crear la H_o para tener una base sobre la que medir una alternativa a nuestra H_a con los datos obervados, calcular un test estadístico desde los datos dados, y finalmente comparar este test estadístico con la distribución hipotética. Esta comparación nos dice cómo de alejado el valor del test está de la hipótesis alternativa. Pvalue es la probabilidad sobre la H_0 de obtener evidencia tanto o más alejada de test estadístico obtenido de los datos observados en la dirección de la hipótesis alternativa.

Por tanto si pvalue es pequeño, o H_o es cierto y los datos observados tienen eventos raros o H_o es falso. Pvalye es el valor más pequeño de alfa en el cual rechazaríamos la H_0

P-value permite elegir un porcentaje de error Typo I , comparando ese porcenaje de error si P-value es menor rechazar la hipotesis nula. Hay muchos mal entendidos en su interpretación, es por eso que tiene muchos defensores y muchos detractores.

Vamos a calcular, cómo de inusal son nuestros dato estimado asumiento que Ho es cierto

qownnotes-media-PC7296

qownnotes-media-PC7296

Para calcular pvalue en una distribución nomal, existe la función pnorm() con los distintos parámetros de la distribución, uno de ellos, lower.trail = TRUE/FALSE le decimos si queremos calcularlo por un sólo extremo o por los dos extremos.En el ejemplo queríamos ver el valor de alfa = 2, si sólo cogemos un extremo tendremos alfa >2

qownnotes-media-H11100

qownnotes-media-H11100

Si la ditribución es binomial pbinom.

qownnotes-media-wSCjdp

qownnotes-media-wSCjdp

qownnotes-media-s11100 qownnotes-media-n11100

Si la distribución es de poisson ppois()

qownnotes-media-HOxTYi

qownnotes-media-HOxTYi

qownnotes-media-u11100

qownnotes-media-u11100

Proporcionando el valor de pvalue en lugar de un alfa (probabilidad de cometer error tipo I) permites a los revisores de tu trabajo realizar hipótesis con el alpha que ellos elijan. La norma general será que si pvalue es menor que el alfa, rechazas H_o.

Si la hipótesis abarca ambos extremos, lo que se suele hacer es doblar el pvalue.

qownnotes-media-CPmAdQ

qownnotes-media-CPmAdQ

qownnotes-media-ptQNCu

qownnotes-media-ptQNCu

qownnotes-media-dQbvin

qownnotes-media-dQbvin

qownnotes-media-tgvovm

qownnotes-media-tgvovm

qownnotes-media-krqBtz

qownnotes-media-krqBtz

qownnotes-media-QHoCgi

qownnotes-media-QHoCgi

qownnotes-media-sNVxow

qownnotes-media-sNVxow

qownnotes-media-rMQxye

qownnotes-media-rMQxye

qownnotes-media-MKVoiM

qownnotes-media-MKVoiM

qownnotes-media-kCHvKO

qownnotes-media-kCHvKO

TEST

qownnotes-media-xOUDhE

qownnotes-media-xOUDhE

qownnotes-media-uOtGOY

qownnotes-media-uOtGOY

qownnotes-media-nLCkjG

qownnotes-media-nLCkjG

qownnotes-media-SaFnex

qownnotes-media-SaFnex

qownnotes-media-ALBoXw

qownnotes-media-ALBoXw

qownnotes-media-cPFSzZ

qownnotes-media-cPFSzZ

qownnotes-media-hbyEwA

qownnotes-media-hbyEwA

PODER DE CONTRASTE

Poder es la probabilidad de rechazar la hipótesis nula cuando es falsa. El error tipo II es llamada Beta, es la proabilidad de no rechazar H_o cuando en realidad es falsa, por tanto el poder = 1-beta.

Se usa por ejemplo, cuando quieres determinar si el tamaño de la muestra es sufientemente grande para arrojar resultados significativos que no se consideren por azar.

En el ejemplo de eventos nocturnos que causan apnea del sueño, la zona azul que queda a la derecha de la línea vertical es el poder de contraste

qownnotes-media-C16176

qownnotes-media-C16176

El poder crecerá cuanto más grande sea la muestra y contra más alejado esté el velor que intentamos inferir del calculado por H_0, en el dibujo de arriba, cuando mayor sea la media de la curva azu, más se aleja del alfa que es la recta y negral mayor área caerá sobre ella a partir de la línea vertical que representa el 5% de H_0 (alfa o error tipo I). En este gráfico vemos cómo según aumenta el valor de media, el poder aumenta, y de manera mucho más incremental si aumentamos el número de muestras

qownnotes-media-sv4036

qownnotes-media-sv4036

Cuanto más pequeá sea la desviación, más delgada será la campana, por lo que las curvas se solapan menos y su poder será mayor. osea poder de contraste decrecerá según aumente la desviación estándar ya que significa más variabilidad en los datos y más se solapen H_0 y H_A

qownnotes-media-Y16176

qownnotes-media-Y16176

qownnotes-media-zvFTEF

qownnotes-media-zvFTEF

Podemos calcular el poder de dos maneras, o bien manual o bien con la funcion t.test.

qownnotes-media-izEfWo

qownnotes-media-izEfWo

Para calcularlo de manera manual, usaremos la función pnorm():

z <- qnorm(.95) pnorm(30+z, mean = 32, lower.tail=FALSE)

Esto nos devolverá 0.64, que signiica que con un 64% podemos rechazar H_o cuando es falsa.

Tenemos la función power.t.test que nos realiza esta tarea para distribuciones Tstudent

qownnotes-media-fJ7448

qownnotes-media-fJ7448

En este caso vemos que tanto Delta (diferencia entre medias) como sd, si se cambian proporcionalmente, el cambio no afecta al poder. Este ration (Delta/sd) se llama effect size.

Se puede utilizar la misma fórmula para saber el número de individuos mínimos necesitamos para cumplir con un power concreto, esto es muy útil para saber el tamaño mínimo de la muestra

qownnotes-media-in7448

qownnotes-media-in7448

EJERCICIOS

Type swirl() when you are ready to begin.

Attaching package: ‘swirl’

The following objects are masked from ‘package:shinyjs’:

info, reset

Warning message: package ‘swirl’ was built under R version 3.3.3

swirl()

Welcome to swirl! Please sign in. If you’ve been here before, use the same name as you did then. If you are new, call yourself something unique.

What shall I call you? miguel

Would you like to continue with one of these lessons?

1: Statistical Inference Hypothesis Testing 2: Statistical Inference Probability2 3: Statistical Inference T Confidence Intervals 4: Statistical Inference Variance 5: No. Let me start something new.

Selection: 5

Please choose a course, or type 0 to exit swirl.

1: Exploratory Data Analysis 2: Getting and Cleaning Data 3: R Programming 4: Statistical Inference 5: Take me to the swirl course repository!

Selection: 4

Please choose a lesson, or type 0 to return to course menu.

1: Introduction 2: Probability1 3: Probability2 4: ConditionalProbability 5: Expectations
6: Variance 7: CommonDistros 8: Asymptotics 9: T Confidence Intervals 10: Hypothesis Testing
11: P Values 12: Power 13: Multiple Testing 14: Resampling

Selection: 12

Attempting to load lesson dependencies…
Package ‘reshape2’ loaded correctly!
Package ‘ggplot2’ loaded correctly!

| | 0%

Power. (Slides for this and other Data Science courses may be found at github https://github.com/DataScienceSpecialization/courses/. If you care to use them, they must be downloaded as a zip file and viewed locally. This lesson corresponds to 06_Statistical_Inference/11_Power.)

|= | 1%

In this lesson, as the name suggests, we’ll discuss POWER, which is the probability of rejecting the null hypothesis when it is false, which is good and proper.

|=== | 2%

Hence you want more POWER.

|==== | 3%

Power comes into play when you’re designing an experiment, and in particular, if you’re trying to determine if a null result (failing to reject a null hypothesis) is meaningful. For instance, you might have to determine if your sample size was big enough to yield a meaningful, rather than random, result.

|====== | 4%

Power gives you the opportunity to detect if your ALTERNATIVE hypothesis is true.

|======= | 5%

Do you recall the definition of a Type II error? Remember, errors are bad.

1: Rejecting a true null hypothesis 2: Miscalculating a t score 3: Misspelling the word hypothesis 4: Accepting a false null hypothesis

Selection: 4

All that practice is paying off!

|========= | 7%

Beta is the probability of a Type II error, accepting a false null hypothesis; the complement of this is obviously (1 - beta) which represents the probability of rejecting a false null hypothesis. This is good and this is POWER!

|========== | 8%

Recall our previous example involving the Respiratory Distress Index and sleep disturbances. Our null hypothesis H_0 was that mu = 30 and our alternative hypothesis H_a was that mu > 30.

|============ | 9%

Which of the following expressions represents our test statistic under this null hypothesis? Here X’ represents the sample mean, s is the sample std deviation, and n is the sample size. Assume X’ follows a t distribution.

1: (X’-30)/(s^2/n) 2: 30/(s/sqrt(n)) 3: (X’-30)/(s/sqrt(n)) 4: X’/(s^2/n)

Selection: 3

Nice work!

|============= | 10%

In the expression for the test statistic (X’-30)/(s/sqrt(n)) what does (s/sqrt(n)) represent?

1: a standard measure 2: a standard error 3: a standard sample 4: a standard variance 5: a standard bearer

Selection: 2

That’s the answer I was looking for.

|=============== | 11%

Suppose we’re testing a null hypothesis H_0 with an alpha level of .05. Since H_a proposes that mu > 30 (the mean hypothesized by H_0), power is the probability that the true mean mu is greater than the (1-alpha) quantile or qnorm(.95). For simplicity, assume we’re working with normal distributions of which we know the variances.

|================ | 12%

Here’s the picture we’ve used a lot in these lessons. As you know, the shaded portion represents 5% of the area under the curve. If a test statistic fell in this shaded portion we would reject H_0 because the sample mean is too far from the mean (center) of the distribution hypothesized by H_0. Instead we would favor H_a, that mu > 30. This happens with probability .05.

|================== | 13%

You might well ask, “What does this have to do with POWER?” Good question. We’ll look at some pictures to show you.

|=================== | 14%

First we have to emphasize a key point. The two hypotheses, H_0 and H_a, actually represent two distributions since they’re talking about means or centers of distributions. H_0 says that the mean is mu_0 (30 in our example) and H_a says that the mean is mu_a.

|===================== | 15%

We’re assuming normality and equal variance, say sigma^2/n, for both hypotheses, so under H_0, X’~ N(mu_0, sigma^2/n) and under H_a, X’~ N(mu_a, sigma^2/n).

|====================== | 16%

Here’s a picture with the two distributions. We’ve drawn a vertical line at our favorite spot, at the 95th percentile of the red distribution. To the right of the line lies 5% of the red distribution.

|======================== | 17%

Quick quiz! Which distribution represents H_0?

1: the blue 2: the red

Selection: 1

Nice try, but that’s not exactly what I was hoping for. Try again.
The two distributions have the same spread (variance). They differ in their means (centers). Which has a mean equal to that hypothesized by H_0?

1: the blue 2: the red

Selection: 2

Keep working like that and you’ll get there!

|========================= | 18%

Which distribution represents H_a?

1: the red 2: the blue

Selection: 2

Excellent work!

|=========================== | 20%

From the picture, what is the mean proposed by H_a?

1: 28 2: 32 3: 36 4: 30

Selection: 2

Nice work!

|============================ | 21%

See how much of the blue distribution lies to the right of that big vertical line?

|============================== | 22%

That, my friend, is POWER!

|=============================== | 23%

It’s the area under the blue curve (H_a) to the right of the vertical line.

|================================= | 24%

Note that the placement of the vertical line depends on the null distribution. Here’s another picture with fatter distributions. The vertical line is still at the 95th percentile of the null (red) distribution and 5% of the distribution still lies to its right. The line is calibrated to mu_0 and the variance.

|================================== | 25%

Back to our original picture.

|==================================== | 26%

We’ve shamelessly stolen plotting code from the slides so you can see H_a in action. Let’s look at pictures before we delve into numbers. We’ve fixed mu_0 at 30, sigma (standard deviation) at 4 and n (sample size) at 16. The function myplot just needs an alternative mean, mu_a, as argument. Run myplot now with an argument of 34 to see what it does.

myplot(34)

All that practice is paying off!

|===================================== | 27%

The distribution represented by H_a moved to the right, so almost all (100%) of the blue curve is to the right of the vertical line, indicating that with mu_a=34, the test is more powerful, i.e., there’s a higher probability that it’s correct to reject the null hypothesis since it appears false. Now try myplot with an argument of 33.3.

myplot(33.3)

You got it right!

|======================================= | 28%

This isn’t as powerful as the test with mu_a=34 but it makes a pretty picture. Now try myplot with an argument of 30.

myplot(30)

That’s the answer I was looking for.

|======================================== | 29%

Uh Oh! Did the red curve disappear? No. it’s just under the blue curve. The power now, the area under the blue curve to the right of the line, is exactly 5% or alpha!

|========================================== | 30%

So what did we learn?

|=========================================== | 32%

First, power is a function that depends on a specific value of an alternative mean, mu_a, which is any value greater than mu_0, the mean hypothesized by H_0. (Recall that H_a specified mu>30.)

|============================================= | 33%

Second, if mu_a is much bigger than mu_0=30 then the power (probability) is bigger than if mu_a is close to 30. As mu_a approaches 30, the mean under H_0, the power approaches alpha.

|============================================== | 34%

Just for fun try myplot with an argument of 28.

myplot(28)

You nailed it! Good job!

|================================================ | 35%

We see that the blue curve has moved to the left of the red, so the area under it, to the right of the line, is less than the 5% under the red curve. This then is even less powerful and contradicts H_a so it’s not worth looking at.

|================================================= | 36%

Here’s a picture of the power curves for different sample sizes. Again, this uses code “borrowed” from the slides. The alternative means, the mu_a’s, are plotted along the horizontal axis and power along the vertical.

|=================================================== | 37%

What does the graph show us about mu_a?

1: power is independent of mu_a 2: as it gets bigger, it gets less powerful 3: as it gets bigger, it gets more powerful

Selection: 3

All that hard work is paying off!

|==================================================== | 38%

What does the graph show us about sample size?

1: as it gets bigger, it gets more powerful 2: as it gets bigger, it gets less powerful 3: power is independent of sample size

Selection: 1

Excellent work!

|====================================================== | 39%

Now back to numbers. Our test for determining rejection of H_0 involved comparing a test statistic, namely Z=(X’-30)/(sigma/sqrt(n)), against some quantile, say Z_95, which depended on our level size alpha (.05 in this case). H_a proposed that mu > mu_0, so we tested if Z>Z_95. This is equivalent to X’ > Z_95 * (sigma/sqrt(n)) + 30, right?

|======================================================= | 40%

Recall that nifty R function pnorm, which gives us the probability that a value drawn from a normal distribution is greater or less than/equal to a specified quantile argument depending on the flag lower.tail. The function also takes a mean and standard deviation as arguments.

|========================================================= | 41%

Suppose we call pnorm with the quantile 30 + Z_95 * (sigma/sqrt(n)) and specify mu_a as our mean argument. This would return a probability which we can interpret as POWER. Why?

|========================================================== | 42%

Recall our picture of two distributions. 30 + Z_95 * (sigma/sqrt(n)) represents the point at which our vertical line falls. It’s the point on the null distribution at the (1-alpha) level.

|============================================================ | 43%

Study this picture. Calling pnorm with 30 + Z_95 * (sigma/sqrt(n)) as the quantile and mu_a, say 32, as the mean and lower.tail=FALSE does what?

1: returns the area under the red curve to the left of the line 2: returns the area under the red curve to the right of the line 3: returns the area under the blue curve to the left of the line 4: returns the area under the blue curve to the right of the line

Selection: 4

Your dedication is inspiring!

|============================================================= | 45%

Let’s try some examples now. Before we do, what do we know pnorm will return if we specify a quantile less than the mean?

1: an answer dependent on beta 2: an answer less than .50 3: an answer dependent on alpha 4: an answer greater than 1

Selection: 2

Excellent work!

|=============================================================== | 46%

First, define a variable z as qnorm(.95)

z <- qnorm(.95)

You are amazing!

|================================================================ | 47%

Run pnorm now with the quantile 30+z, mean=30, and lower.tail=FALSE. We’ve specified sigma and n so that the standard deviation of the sample mean is 1.

z <- pnorm(30+z,mean=30,lower.tail = FALSE)

Not quite, but you’re learning! Try again. Or, type info() for more options.
Type pnorm(30+z,mean=30,lower.tail=FALSE) at the command prompt.

pnorm(30+z,mean=30,lower.tail = FALSE) [1] 0.05

Your dedication is inspiring!

|================================================================== | 48%

That’s not surprising, is it? With the mean set to mu_0 the two distributions, null and alternative, are the same and power=alpha. Now run pnorm now with the quantile 30+z, mean=32, and lower.tail=FALSE.

pnorm(30+z,mean=32,lower.tail = FALSE) [1] 0.63876

You are really on a roll!

|=================================================================== | 49%

See how this is much more powerful? 64% as opposed to 5%. When the sample mean is quite different from (many standard errors greater than) the mean hypothesized by the null hypothesis, the probability of rejecting H_0 when it is false is much higher. That is power!

|==================================================================== | 50%

Let’s look again at the portly distributions.

|====================================================================== | 51%

With this standard deviation=2 (fatter distribution) will power be greater or less than with the standard deviation=1?

1: less than 2: the same 3: greater

Selection: 1

Keep up the great work!

|======================================================================= | 52%

To see this, run pnorm now with the quantile 30+z, mean=32 and sd=1. Don’t forget to set lower.tail=FALSE so you get the right tail.

pnorm(30+z,mean=32,sd=1,lower.tail = FALSE) [1] 0.63876

Excellent job!

|========================================================================= | 53%

Now run pnorm now with the quantile 30+z*2, mean=32 and sd=2. Don’t forget to set lower.tail=FALSE so you get the right tail.

pnorm(30+z,mean=32,sd=2,lower.tail = FALSE) [1] 0.5704709

Not exactly. Give it another go. Or, type info() for more options.
Type pnorm(30+z*2,mean=32,sd=2,lower.tail=FALSE) at the command prompt.

pnorm(30+z*2,mean=32,sd=2,lower.tail = FALSE) [1] 0.259511

Excellent job!

|========================================================================== | 54%

See the power drain from 64% to 26% ? Let’s review some basic facts about power. We saw before in our pictures that the power of the test depends on mu_a. When H_a specifies that mu > mu_0, then as mu_a grows and exceeds mu_0 increasingly, what happens to power?

1: it doesn’t change 2: it decreases 3: it increases

Selection: 3

You are quite good my friend!

|============================================================================ | 55%

Here’s another question. Recall our power curves from before.

|============================================================================= | 57%

As the sample size increases, what happens to power?

1: it decreases 2: it increases 3: it doesn’t change

Selection: 2

You nailed it! Good job!

|=============================================================================== | 58%

Here’s another one. More power curves.

|================================================================================ | 59%

As variance increases, what happens to power?

1: it doesn’t change 2: it increases 3: it decreases

Selection: 3

All that hard work is paying off!

|================================================================================== | 60%

Here’s another one. And even more power curves.

|=================================================================================== | 61%

As alpha increases, what happens to power?

1: it decreases 2: it doesn’t change 3: it increases

Selection: 3

Keep up the great work!

|===================================================================================== | 62%

If H_a proposed that mu != mu_0 we would calculate the one sided power using alpha / 2 in the direction of mu_a (either less than or greater than mu_0). (This is only approximately right, it excludes the probability of getting a large test statistic in the opposite direction of the truth.

|====================================================================================== | 63%

Since power goes up as alpha gets larger would the power of a one-sided test be greater or less than the power of the associated two sided test?

1: less than 2: greater 3: they’re the same

Selection: 3

Not quite! Try again.
The quantity alpha is bigger than alpha/2 so it’s got more power.

1: less than 2: greater 3: they’re the same

Selection: 2

All that hard work is paying off!

|======================================================================================== | 64%

Finally, if H_a specified that mu < mu_0 could we still do the same kind of power calculations?

1: No 2: Yes

Selection: 2

Nice work!

|========================================================================================= | 65%

Suppose H_a says that mu > mu_0. Then power = 1 - beta = Prob ( X’ > mu_0 + z_(1-alpha) * sigma/sqrt(n)) assuming that X’~ N(mu_a,sigma^2/n). Which quantities do we know in this statement, given the context of the problem? Let’s work through this.

|=========================================================================================== | 66%

What does the null hypothesis H_0 tell us that the population mean equals?

1: alpha 2: mu_0 3: beta 4: mu_a

Selection: 2

All that hard work is paying off!

|============================================================================================ | 67%

After the null mean mu_0 is proposed what does the designer of the hypothesis test specify in order to reject or fail-to-reject H_0? In other words, what is the level size of the test?

1: mu_0 2: alpha 3: mu_a 4: beta

Selection: 2

You nailed it! Good job!

|============================================================================================== | 68%

So we know that the quantities mu_0 and alpha are specified by the test designer. In the statement 1 - beta = Prob( X’ > mu_0 + z_(1-alpha) * sigma/sqrt(n)) given mu_a > mu_0, mu_0 and alpha are specified, and X’ depends on the data. The other four quantities, (beta, sigma, n, and mu_a), are all unknown.

|=============================================================================================== | 70%

It should be obvious that specifying any three of these unknowns will allow us to solve for the missing fourth. Usually, you only try to solve for power (1-beta) or the sample size n.

|================================================================================================= | 71%

An interesting point is that power doesn’t need mu_a, sigma and n individually. Instead only sqrt(n)*(mu_a - mu_0) /sigma is needed. The quantity (mu_a - mu_0) / sigma is called the EFFECT SIZE. This is the difference in the means in standard deviation units. It is unit free so it can be interpreted in different settings.

|================================================================================================== | 72%

We’ll work through some examples of this now. However, instead of assuming that we’re working with normal distributions let’s work with t distributions. Remember, they’re pretty close to normal with large enough sample sizes.

|==================================================================================================== | 73%

Power is still a probability, namely P( (X’ - mu_0)/(S /sqrt(n)) > t_(1-alpha, n-1) given H_a that mu > mu_a ). Notice we use the t quantile instead of the z. Also, since the proposed distribution is not centered at mu_0, we have to use the non-central t distribution.

|===================================================================================================== | 74%

R comes to the rescue again with the function power.t.test. We can omit one of the arguments and the function solves for it. Let’s first use it to solve for power.

|======================================================================================================= | 75%

We’ll run it three times with the same values for n (16) and alpha (.05) but different delta and standard deviation values. We’ll show that if delta (difference in means) divided by the standard deviation is the same, the power returned will also be the same. In other words, the effect size is constant for all three of our tests.

|======================================================================================================== | 76%

We’ll specify a positive delta; this tells power.t.test that H_a proposes that mu > mu_0 and so we’ll need a one-sided test. First run power.t.test(n = 16, delta = 2 / 4, sd=1, type = “one.sample”, alt = “one.sided”)$power .

power.t.test(n = 16, delta = 2 / 4, sd=1, type = “one.sample”,alt = “one.sided”)$power [1] 0.6040329

You are doing so well!

|========================================================================================================== | 77%

Now change delta to 2 and sd to 4. Keep everything else the same.

power.t.test(n = 16, delta = 2, sd=4, type = “one.sample”,alt = “one.sided”)$power [1] 0.6040329

Excellent job!

|=========================================================================================================== | 78%

Same answer, right? Now change delta to 100 and sd to 200. Keep everything else the same.

power.t.test(n = 16, delta = 100, sd=200, type = “one.sample”,alt = “one.sided”)$power [1] 0.6040329

Nice work!

|============================================================================================================= | 79%

So keeping the effect size (the ratio delta/sd) constant preserved the power. Let’s try a similar experiment except now we’ll specify a power we want and solve for the sample size n.

|============================================================================================================== | 80%

First run power.t.test(power = .8, delta = 2 / 4, sd=1, type = “one.sample”, alt = “one.sided”)$n .

power.t.test(power = .8, delta = 2 / 4, sd=1, type = “one.sample”, alt = “one.sided”)$n [1] 26.13751

Nice work!

|================================================================================================================ | 82%

Now change delta to 2 and sd to 4. Keep everything else the same.

power.t.test(power = .8, delta = 2, sd=4, type = “one.sample”, alt = “one.sided”)$n [1] 26.13751

Keep up the great work!

|================================================================================================================= | 83%

Same answer, right? Now change delta to 100 and sd to 200. Keep everything else the same.

power.t.test(power = .8, delta = 100, sd=200, type = “one.sample”, alt = “one.sided”)$n [1] 26.13751

Excellent job!

|=================================================================================================================== | 84%

Now use power.t.test to find delta for a power=.8 and n=26 and sd=1

power.t.test(power = .8, n = 26, sd=1, type = “one.sample”, alt = “one.sided”)

 One-sample t test power calculation

          n = 26
      delta = 0.5013986
         sd = 1
  sig.level = 0.05
      power = 0.8
alternative = one.sided
You’re close…I can feel it! Try it again. Or, type info() for more options.
Type power.t.test(power = .8, n=26, sd=1, type = “one.sample”, alt = “one.sided”)$delta.

power.t.test(power = .8, n = 26, sd=1, type = “one.sample”, alt = “one.sided”)$delta [1] 0.5013986

Great job!

|==================================================================================================================== | 85%

Not a surprising result, is it? It told you before that with an effect size of .5 and power .8, you need a sample size a little more than 26. Now run it with n=27.

power.t.test(power = .8, n = 27, sd=1, type = “one.sample”, alt = “one.sided”)$delta [1] 0.4914855

You’re the best!

|====================================================================================================================== | 86%

What do you think will happen if you doubled sd to 2 and ran the same test?

1: delta will halve 2: delta will double 3: delta won’t change

Selection: 1

Keep trying!
Since you’re doubling the denominator (sd) you have to double the numerator (delta) in order to keep the effect size constant.

1: delta won’t change 2: delta will halve 3: delta will double

Selection: 1

Not quite! Try again.
Since you’re doubling the denominator (sd) you have to double the numerator (delta) in order to keep the effect size constant.

1: delta will halve 2: delta won’t change 3: delta will double

Selection: 3

Great job!

|======================================================================================================================= | 87%

Now for a quick review. We call this the power.u.test since it comes after the power.t.test. LOL.

|========================================================================================================================= | 88%

1. The level of a test is specified by what?

1: gamma 2: None of the others 3: alpha 4: beta 5: delta

Selection: 5

One more time. You can do it!
The level refers to the probability of rejecting the null hypothesis when it is true. This is the first Greek letter.

1: gamma 2: delta 3: None of the others 4: alpha 5: beta

Selection: 4

Perseverance, that’s the answer.

|========================================================================================================================== | 89%

2. What is a Type II error?

1: accepting a true hypothesis 2: rejecting a false hypothesis 3: rejecting a true hypothesis 4: accepting a false hypothesis

Selection: 4

All that hard work is paying off!

|============================================================================================================================ | 90%

3. What is power?

1: delta 2: thrilling 3: None of the others 4: alpha 5: beta 6: gamma

Selection: 3

Perseverance, that’s the answer.

|============================================================================================================================= | 91%

4. You’re a perfectionist designing an experiment and you want both alpha and beta to be small. Can they both be 0 for this single test?

1: No 2: Yes

Selection: 1

Excellent job!

|=============================================================================================================================== | 92%

5. Suppose H_0 proposes mu = mu_0 and H_a proposes that mu < mu_0. You’ll test a series of mu_a with power != alpha. Which of the following is NOT true?

1: huh? 2: mu_a-mu_0 < 0 3: mu_a-mu_0=0 4: mu_0-mu_a > 0

Selection: 3

Keep up the great work!

|================================================================================================================================ | 93%

6. Suppose H_0 proposes mu = mu_0 and H_a proposes that mu < mu_0. Which of the following is true?

1: mu_0=mu_a maximizes the power 2: the smaller mu_0-mu_a the more powerful the test 3: the smaller mu_a-mu_0 the more powerful the test

Selection: 2

That’s not exactly what I’m looking for. Try again.
Here mu_a < mu_0 and the smaller mu_a-mu_0 is, the easier it is to discriminate between mu_a and mu_0.

1: the smaller mu_0-mu_a the more powerful the test 2: mu_0=mu_a maximizes the power 3: the smaller mu_a-mu_0 the more powerful the test

Selection: 3

You are doing so well!

|================================================================================================================================== | 95%

7. Which expression represents the size effect?

1: (mu_a - mu_0) / sqrt(sigma) 2: (mu_a - mu_0) / sqrt(n) 3: (mu_a - mu_0) / n 4: (mu_a - mu_0) / sigma

Selection: 4

You are doing so well!

|=================================================================================================================================== | 96%

8. True or False? More power is better than less power.

1: False 2: True

Selection: 2

That’s correct!

|===================================================================================================================================== | 97%

9. True or False? A larger beta (call it beta_max) is more powerful than a smaller beta.

1: True 2: False

Selection: 2

You’re the best!

|====================================================================================================================================== | 98%

10. True or False? The larger the sample size the less powerful the test.

1: False 2: True

Selection: 1

That’s correct!

|======================================================================================================================================== | 99%

Congrats! You finished this powerful lesson. We hope you feel emPOWERED.

|=========================================================================================================================================| 100%

Would you like to receive credit for completing this course on Coursera.org?

1: No 2: Yes

COMPARACIONES MÚLTIPLES (MULTIPLE TESTING)

Es un todo un subcampo detro de la inferencia estadística. Trata de reducir y ajustar los errores tipo I y II para evitar falsos positivos (rechazar H_o cuando no se debe) y falsos negativos (rechazar H_a cuando no se debe).

Imaginemos que queremos probar una hipótesis sobre múltiples grupos. Lo que haremos será comparar las medias obtenidas de todos esos grupos.

Un ejemplo claro que podemos tener en esta situación es que, pongamos que hacemos 10.000 test, si calculamos un PValue de manera que P< alfa estaremos controlando los falsos positivos en un nivel de media alfa. Suponemos P< 0.05 significante, por lo que 10000*0.05 = 500 posibles falsos positivos, que parece un número muy alto. Cómo evitamos esto?

Vamos a jugar con ALFA equilibrando la reducción del número de falsos positivos con que tengamos una ALFA muy baja y no podamos rechazar nunca Ho

qownnotes-media-y17684

qownnotes-media-y17684

FWER Para controlar el falso positivo: Family wise error rate (FWER) que es la probabilidad de encontrar al menos 1 falso positivo. Se puede intentar controlar este FWER con la corrección de Bonferroni. Lo bueno de este método es que es muy rápido, pero baja tanto ALFA que nos será muy dificil rechazar Ho

qownnotes-media-X17684

qownnotes-media-X17684

FDR

Otro método para controlar los falsos positivos es False discovery rate (FDR) most popular in genomics E(V/R) = Ratio de falsos positivos (número de falsos positivos divido número de positivos). En este caso se aplica la corrección BH (Benjamini-Hochberg). Es fácil de calcular y aunque permite más falsos positivos es menos conservador y permite rechazar más facilemnte Ho

qownnotes-media-Z17684

qownnotes-media-Z17684

En esta gráfica, se ve que el número de resultados significantes sin corrección son 4, con la corrección HC baja a 3 y con la de bonferoni a 2

qownnotes-media-F17684

qownnotes-media-F17684

Adjusted Pvalues

Hay otra tercera aproximación que es en lugar de buscar un umbral a ALFA ajustar pvalues de las muestras

qownnotes-media-Y17684

qownnotes-media-Y17684

qownnotes-media-fMlJWx

qownnotes-media-fMlJWx

qownnotes-media-KNlCiF

qownnotes-media-KNlCiF

qownnotes-media-WpNVqa

qownnotes-media-WpNVqa

qownnotes-media-cGLCMd

qownnotes-media-cGLCMd

qownnotes-media-ubeOBL

qownnotes-media-ubeOBL

qownnotes-media-RrJxGj

qownnotes-media-RrJxGj

En R hay una función llamada p.adjust que hace esto, teniendo un vector de mil pvalues:

qownnotes-media-d14368 qownnotes-media-d14368 qownnotes-media-c14368 qownnotes-media-r14368 qownnotes-media-w14368 qownnotes-media-Z14368 qownnotes-media-n14368 qownnotes-media-R14368 qownnotes-media-c14368

EJERCICIOS

Given that data is valuable and we’d like to get the most out of it, we might use it to test several hypotheses. If we have an alpha level of .05 and we test 20 hypotheses, then on average, we expect one error, just by chance.

|===== | 5%

Another potential problem is that after running several tests, only the lowest p-value might be reported OR all p-values under some threshold might be considered significant. Undoubtedly, some of these would be false.

|======= | 7%

Luckily, we have clever ways of minimizing errors in this situation. That’s what we’ll address. We’ll define specific error measures and then statistical ways of correcting or limiting them.

|========= | 8%

Multiple testing is particularly relevant now in this age of BIG data. Statisticians are tasked with questions such as “Which variables matter among the thousands measured?” and “How do you relate unrelated information?”

|========== | 10%

Since multiple testing addresses compensating for errors let’s review what we know about them. A Type I error is

1: rejecting a false hypothesis 2: failing to reject a true hypothesis 3: rejecting a true hypothesis 4: failing to reject a false hypothesis

Selection: 4

Nice try, but that’s not exactly what I was hoping for. Try again.
Eliminate the two choices that are not errors. A Type I error involves rejection.

1: rejecting a true hypothesis 2: failing to reject a true hypothesis 3: rejecting a false hypothesis 4: failing to reject a false hypothesis

Selection: 3

Keep trying!
Eliminate the two choices that are not errors. A Type I error involves rejection.

1: rejecting a true hypothesis 2: failing to reject a false hypothesis 3: failing to reject a true hypothesis 4: rejecting a false hypothesis

Selection: 1

All that practice is paying off!

|============ | 11%

In an American court, an example of a Type I error is

1: convicting an innocent person 2: letting the indicted off on a technicality 3: acquitting a guilty person

Selection: 3

Keep trying!
In an American court, the null hypothesis is that the accused is innocent. If he is convicted when he really is innocent then the null hypothesis is rejected incorrectly.

1: letting the indicted off on a technicality 2: acquitting a guilty person 3: convicting an innocent person

Selection: 3

You’re the best!

|============== | 13%

A Type II error is

1: rejecting a false hypothesis 2: failing to reject a false hypothesis 3: failing to reject a true hypothesis 4: rejecting a true hypothesis

Selection: 1

Almost! Try again.
Eliminate the two choices that are not errors. A Type II error involves failing to reject.

1: failing to reject a false hypothesis 2: rejecting a true hypothesis 3: rejecting a false hypothesis 4: failing to reject a true hypothesis

Selection: 1

That’s correct!

|=============== | 15%

In an American court, an example of a Type II error is

1: letting the indicted off on a technicality 2: convicting an innocent person 3: acquitting a guilty person

Selection: 3

Great job!

|================= | 16%

Good. Let’s continue reviewing. The null hypothesis

1: tells us the origins of the number 0 2: represents the status_quo and is assumed true 3: is a big nothing that statisticians like to gossip about 4: is never true

Selection: 2

That’s the answer I was looking for.

|=================== | 18%

The p-value is “the probability under the null hypothesis of obtaining evidence as or more extreme than your test statistic (obtained from your observed data) in the direction of the alternative hypothesis." Of course p-values are related to significance or alpha levels, which are set before the test is conducted (often at 0.05).

|==================== | 20%

If a p-value is found to be less than alpha (say 0.05), then the test result is considered statistically significant, i.e., surprising and unusual, and the null hypothesis (the status quo) is ?

1: revised 2: renamed the aleph null hypothesis 3: rejected 4: accepted

Selection: 3

You got it right!

|====================== | 21%

Now consider this chart copied from http://en.wikipedia.org/wiki/Familywise_error_rate. Suppose we’ve tested m null hypotheses, m_0 of which are actually true, and m-m_0 are actually false. Out of the m tests R have been declared significant, that is, the associated p-values were less than alpha, and m-R were nonsignificant, or boring results.

|======================== | 23%

Looking at the chart, which variables are known?

1: A,B,C 2: S,T,U,V 3: m_0, and m 4: m and R

Selection: 4

That’s a job well done!

|========================== | 25%

In testing the m_0 true null hypotheses, V results were declared significant, that is, these tests favored the alternative hypothesis. What type of error does this represent?

1: Type I 2: Type III 3: Type II 4: a serious one

Selection: 1

Excellent work!

|=========================== | 26%

Another name for a Type I error is False Positive, since it is falsely claiming a significant (positive) result.

|============================= | 28%

Of the m-m_0 false null hypotheses, T were declared nonsignificant. This means that these T null hypotheses were accepted (failed to be rejected). What type of error does this represent?

1: Type III 2: Type II 3: a serious one 4: Type I

Selection: 2

Perseverance, that’s the answer.

|=============================== | 30%

Another name for a Type II error is False Negative, since it is falsely claiming a nonsignificant (negative) result.

|================================ | 31%

A rose by any other name, right? Consider the fraction V/R.

|================================== | 33%

The observed R represents the number of test results declared significant. These are ‘discoveries’, something different from the status quo. V is the number of those falsely declared significant, so V/R is the ratio of FALSE discoveries. Since V is a random variable (i.e., unknown until we do an experiment) we call the expected value of the ratio, E(V/R), the False Discovery Rate (FDR).

|==================================== | 34%

A rose by any other name, right? How about the fraction V/m_0? From the chart, m_0 represents the number of true H_0’s and m_0 is unknown. V is the number of those falsely declared significant, so V/m_0 is the ratio of FALSE positives. Since V is a random variable (i.e., unknown until we do an experiment) we call the expected value of the ratio, E(V/m_0), the FALSE POSITIVE rate.

|====================================== | 36%

Another good name for the false positive rate would be

1: false alarm rate 2: the Type II rate 3: a thorn 4: a rose

Selection: 1

You are really on a roll!

|======================================= | 38%

The false positive rate would be closely related to

1: the Type I error rate 2: the Type II error rate 3: a thorny rose

Selection: 1

Perseverance, that’s the answer.

|========================================= | 39%

We call the probability of at least one false positive, Pr(V >= 1) the Family Wise Error Rate (FWER).

|=========================================== | 41%

So how do we control the False Positive Rate?

|============================================ | 43%

Suppose we’re really smart, calculate our p-values correctly, and declare all tests with p < alpha as significant. This means that our false positive rate is at most alpha, on average.

|============================================== | 44%

Suppose we perform 10,000 tests and alpha = .05. How many false positives do we expect on average?

1: 50000 2: 5000 3: 50 4: 500

Selection: 4

Keep working like that and you’ll get there!

|================================================ | 46%

You got it! 500 false positives seems like a lot. How do we avoid so many?

|================================================= | 48%

We can try to control the family-wise error rate (FWER), the probability of at least one false positive, with the Bonferroni correction, the oldest multiple testing correction.

|=================================================== | 49%

It’s very straightforward. We do m tests and want to control the FWER at level alpha so that Pr(V >= 1) < alpha. We simply reduce alpha dramatically. Set alpha_fwer to be alpha/m. We’ll only call a test result significant if its p-value < alpha_fwer.

|===================================================== | 51%

Sounds good, right? Easy to calculate. What would be a drawback with this method?

1: too many results will fail 2: too many results will pass 3: requires too much math

Selection: 1

That’s the answer I was looking for.

|======================================================= | 52%

Another way to limit the false positive rate is to control the false discovery rate (FDR). Recall this is E(V/R). This is the most popular correction when performing lots of tests. It’s used in lots of areas such as genomics, imaging, astronomy, and other signal-processing disciplines.

|======================================================== | 54%

Again, we’ll do m tests but now we’ll set the FDR, or E(V/R) at level alpha. We’ll calculate the p-values as usual and order them from smallest to largest, p_1, p_2,…p_m. We’ll call significant any result with p_i <= (alpha*i)/m. This is the Benjamini-Hochberg method (BH). A p-value is compared to a value that depends on its ranking.

|========================================================== | 56%

This is equivalent to finding the largest k such that p_k <= (k * alpha)/m, (for a given alpha) and then rejecting all the null hypotheses for i=1,…,k.

|============================================================ | 57%

Like the Bonferroni correction, this is easy to calculate and it’s much less conservative. It might let more false positives through and it may behave strangely if the tests aren’t independent.

|============================================================= | 59%

Now consider this chart copied from the slides. It shows the p-values for 10 tests performed at the alpha=.2 level and three cutoff lines. The p-values are shown in order from left to right along the x-axis. The red line is the threshold for No Corrections (p-values are compared to alpha=.2), the blue line is the Bonferroni threshold, alpha=.2/10 = .02, and the gray line shows the BH correction. Note that it is not horizontal but has a positive slope as we expect.

|=============================================================== | 61%

With no correction, how many results are declared significant?

1: 8 2: 6 3: 2 4: 4

Selection: 4

Excellent job!

|================================================================= | 62%

With the Bonferroni correction, how many tests are declared significant?

1: 6 2: 2 3: 4 4: 8

Selection: 2

All that practice is paying off!

|================================================================== | 64%

So the Bonferroni passed only half the results that the No Correction (comparing p-values to alpha) method passed. Now look at the BH correction. How many tests are significant with this scale?

1: 1 2: 7 3: 5 4: 3

Selection: 4

You’re the best!

|==================================================================== | 66%

So the BH correction which limits the FWER is between the No Correction and the Bonferroni. It’s more conservative (fewer significant results) than the No Correction but less conservative (more significant results) than the Bonferroni. Note that with this method the threshold is proportional to the ranking of the values so it slopes positively while the other two thresholds are flat.

|====================================================================== | 67%

Notice how both the Bonferroni and BH methods adjusted the threshold (alpha) level of rejecting the null hypotheses. Another equivalent corrective approach is to adjust the p-values, so they’re not classical p-values anymore, but they can be compared directly to the original alpha.

|======================================================================== | 69%

Suppose the p-values are p_1, … , p_m. With the Bonferroni method you would adjust these by setting p’_i = max(m * p_i, 1) for each p-value. Then if you call all p’_i < alpha significant you will control the FWER.

|========================================================================= | 70%

To demonstrate some of these concepts, we’ve created an array of p-values for you. It is 1000-long and the result of a linear regression performed on random normal x,y pairs so there is no true significant relationship between the x’s and y’s.

|=========================================================================== | 72%

Use the R command head to see the first few entries of the array pValues.

head(pValues) [1] 0.5334915 0.2765785 0.8380943 0.6721730 0.8122037 0.4078675

Great job!

|============================================================================= | 74%

Now count the number of entries in the array that are less than the value .05. Use the R command sum, and the appropriate Boolean expression.

sum() [1] 0

That’s not the answer I was looking for, but try again. Or, type info() for more options.
Type sum(pValues < 0.05) at the command prompt.

sum(pValues < 0.05) [1] 51

You are doing so well!

|============================================================================== | 75%

So we got around 50 false positives, just as we expected (.05*1000=50). The beauty of R is that it provides a lot of built-in statistical functionality. The function p.adjust is one example. The first argument is the array of pValues. Another argument is the method of adjustment. Once again, use the R function sum and a boolean expression using p.adjust with method=“bonferroni” to control the FWER.

p.adjust(pValue,method=“bonferroni”) Error in p.adjust(pValue, method = “bonferroni”) : object ‘pValue’ not found p.adjust(pValues,method=“bonferroni”)

sum(p.adjust(pValues,method=“bonferroni”) < 0.05) [1] 0

Your dedication is inspiring!

|================================================================================ | 77%

So the correction eliminated all the false positives that had passed the uncorrected alpha test. Repeat the same experiment, this time using the method “BH” to control the FDR.

sum(p.adjust(pValues,method=“BH”) < 0.05) [1] 0

Perseverance, that’s the answer.

|================================================================================== | 79%

So the BH method also eliminated all the false positives. Now we’ve generated another 1000-long array of p-values, this one called pValues2. In this data, the first half ( 500 x/y pairs) contains x and y values that are random and the second half contain x and y pairs that are related, so running a linear regression model on the 1000 pairs should find some significant (not random) relationship.

|==================================================================================== | 80%

We also created a 1000-long array of character strings, trueStatus. The first 500 entries are “zero” and the last are “not zero”. Use the R function tail to look at the end of trueStatus.

tail(trueStatus) [1] “not zero” “not zero” “not zero” “not zero” “not zero” “not zero”

Great job!

|===================================================================================== | 82%

Once again we can use R’s greatness to count and tabulate for us. We can call the R function table with two arguments, a boolean such as pValues2<.05, and the array trueStatus. The boolean obviously has two outcomes and each entry of trueStatus has one of two possible values. The function table aligns the two arguments and counts how many of each combination (TRUE,“zero”), (TRUE,“not zero”), (FALSE,“zero”), and (FALSE,“not zero”) appear. Try it now.

table(trueStatus) trueStatus not zero zero 500 500

That’s not the answer I was looking for, but try again. Or, type info() for more options.
Type table(pValues2 < 0.05, trueStatus) at the command prompt.

table(pValues2 < 0.05, trueStatus) trueStatus not zero zero FALSE 0 476 TRUE 500 24

You are doing so well!

|======================================================================================= | 84%

We see that without any correction all 500 of the truly significant (nonrandom) tests were correctly identified in the “not zero” column. In the zero column (the truly random tests), however, 24 results were flagged as significant.

|========================================================================================= | 85%

What is the percentage of false positives in this test?

4.8 [1] 4.8

Give it another try. Or, type info() for more options.
Divide 24 by 500 to get the percentage.

24/500 [1] 0.048

Your dedication is inspiring!

|========================================================================================== | 87%

Just as we expected - around 5% or .05*100.

|============================================================================================ | 89%

Now run the same table function, however, this time use the call to p.adjust with the “bonferroni” method in the boolean expression. This will control the FWER.

table(pValues2 < 0.05, trueStatus) trueStatus not zero zero FALSE 0 476 TRUE 500 24

Nice try, but that’s not exactly what I was hoping for. Try again. Or, type info() for more options.
Type table(p.adjust(pValues2,method=“bonferroni”) < 0.05, trueStatus) at the command prompt.

table(p.adjust(pValues2,method=“bonferroni”) < 0.05, trueStatus) trueStatus not zero zero FALSE 23 500 TRUE 477 0

You nailed it! Good job!

|============================================================================================== | 90%

Since the Bonferroni correction method is more conservative than just comparing p-values to alpha all the truly random tests are correctly identified in the zero column. In other words, we have no false positives. However, the threshold has been adjusted so much that 23 of the truly significant results have been misidentified in the not zero column.

|=============================================================================================== | 92%

Now run the same table function one final time. Use the call to p.adjust with “BH” method in the boolean expression. This will control the false discovery rate.

table(p.adjust(pValues2,method=“BH”) < 0.05, trueStatus) trueStatus not zero zero FALSE 0 487 TRUE 500 13

Great job!

|================================================================================================= | 93%

Again, the results are a compromise between the No Corrections and the Bonferroni. All the significant results were correctly identified in the “not zero” column but in the random (“zero”) column 13 results were incorrectly identified. These are the false positives. This is roughly half the number of errors in the other two runs.

|=================================================================================================== | 95%

Here’s a plot of the two sets of adjusted p-values, Bonferroni on the left and BH on the right. The x-axis indicates the original p-values. For the Bonferroni, (adjusting by multiplying by 1000, the number of tests), only a few of the adjusted values are below 1. For the BH, the adjusted values are slightly larger than the original values.

|===================================================================================================== | 97%

We’ll conclude by saying that multiple testing is an entire subfield of statistical inference. Usually a basic Bonferroni/BH correction is good enough to eliminate false positives, but if there is strong dependence between tests there may be problems. Another correction method to consider is “BY”.

|====================================================================================================== | 98%

Congrats! We hope you liked the multiple concepts and questions you saw in this lesson.

|========================================================================================================| 100%

RESAMPLING INFERENCE

Hasta ahora hemos visto métodos de inferencia estadística para obtener los parámetros de una población en base a su distrubución, pero ¿qué pasa si desconocemos su distribución?¿Cómo podemos obterner la media y varianza poblacional por ejemplo al tirar un dado si sabemos que está trucado? De lo que se trata entonces es de realizar inferencias sobre nuestros datos tomando muestras de nuestros propios datos, es decir haciendo múltiples selecciones aleatorias sobre los datos que hemos recogido, para poder obterner la media de sus medias.

BOOTSTRAPING

Es una técnica que permite simplificar la obtención de los intervalos de confianza en lugar de con métodos matemáticos, calcularlos a través de múltiples muestras aleatorias tomadas de las muestras sin sacar del espacio muestral en cada iteración la muesta obtenida (reemplazo).

Es una técnica no paramétrica y nos serviría para poder detallar en nuestro estudio un intervalo de confianza y el error estándar sobre una única medición tomada, iterando una y otra vez sobre esa muesta para obtener de manera aleatoria uno se sus valores y así simular la población.

qownnotes-media-JkhSIH

qownnotes-media-JkhSIH

qownnotes-media-tanpiX

qownnotes-media-tanpiX

qownnotes-media-oTTWHa

qownnotes-media-oTTWHa

qownnotes-media-HmdAlZ

qownnotes-media-HmdAlZ

qownnotes-media-MvzZKg

qownnotes-media-MvzZKg

qownnotes-media-ugudYv

qownnotes-media-ugudYv

qownnotes-media-CGFVKT

qownnotes-media-CGFVKT

qownnotes-media-zblFDp

qownnotes-media-zblFDp

PERMUTATION TEST

Esta es otra técnica de inferencia por muestreo, en la que en lugar de utilizar T-Student o Chicuadrado, como no sabemos la distribución de los datos comparamos los grupos haciendo permutaciones aleatorias entre ellos mismos y lo comparamos con los resultados más extremos que podemos obtener de su comparativa inicial.

A grandes rasgos, se restan las medias de ambos grupos y se comprueba con otro grupo de permutaciones realizadas cuántas restas de las medias de las permutaciones aleatorias supera la medida inicial, y así saber cuán extremo es la diferencia.

qownnotes-media-WVBDaB

qownnotes-media-WVBDaB

qownnotes-media-lrLfbM

qownnotes-media-lrLfbM

qownnotes-media-RRTRbe

qownnotes-media-RRTRbe

qownnotes-media-BxnNue

qownnotes-media-BxnNue

qownnotes-media-txJbfA

qownnotes-media-txJbfA

EJERCICIOS

The bootstrap is a handy tool for making statistical inferences. It is used in constructing confidence intervals and calculating standard errors for statistics that might be difficult for some reason (e.g., lack of data or no closed form). Wikipedia tells us that bootstrapping is a technique which “allows estimation of the sampling distribution of almost any statistic using very simple methods." Simple is good, right?

|==== | 4%

The beauty of bootstrapping is that it avoids complicated mathematics and instead uses simulation and computation to infer distributional properties you might not otherwise be able to determine.

|====== | 6%

It’s relatively new, developed in 1979, by Bradley Efron, a Stanford statistician. The basic bootstrap principle uses OBSERVED data to construct an ESTIMATED population distribution using random sampling with replacement. From this distribution (constructed from the observed data) we can estimate the distribution of the statistic we’re interested in.

|======= | 7%

So, in bootstrapping the observed data substitutes for what?

1: a statistic 2: a hypothesis 3: a population 4: observations

Selection: 1

Not quite right, but keep trying.
We’ll used the observed data and sample it with replacement, just as we would do to find out information about some population.

1: a statistic 2: a population 3: observations 4: a hypothesis

Selection: 2

You nailed it! Good job!

|========= | 8%

So, in bootstrapping if the observed data is the population, what would the random samplings correspond to?

1: a statistic 2: a hypothesis 3: observations 4: a population

Selection: 3

You got it right!

|========== | 10%

In effect, the original observed sample substitutes for the population. Our samplings become observations from which we estimate a statistic and get an idea about its distribution. This lets us better understand the underlying population (from which we didn’t have enough data).

|============ | 11%

Here’s a critical point. In constructing the estimated distribution we sample the observed data WITH replacement. If the original sample is n long and we sampled n times without replacement what would we get?

1: an entirely new sample 2: a worse sample 3: the original sample permuted 4: a better sample

Selection: 3

Great job!

|============= | 12%

The motivating example from the slides involves computing the average of 50 rolls of a die. Of course we can do this theoretically when we know that the die is fair. Remember, E(x) = Sum(x*p(x)) for x=1,2,…6, and p(x)=1/6 for all values of x.

|============== | 14%

For the heck of it, compute the expected die roll for a fair die.

3 [1] 3

Give it another try. Or, type info() for more options.
Type sum(1:6)/6 at the command prompt.

sum(1:6)/6 [1] 3.5

Perseverance, that’s the answer.

|================ | 15%

Theoretically, the average is 3.5. Here, we’ve run code and plotted a histogram after we took 1000 such averages, each of 50 dice rolls. Note the unusual y-axis scale. We’re displaying this as a density function so the area of the salmon-colored region is theoretically 1. With this scale, though, all the heights of the bins actually add up to 5. So you have to multiply each height by .2 and add up all the results to get 1.

|================= | 17%

The point is, the empirical matches the theoretical. Yay! The highest bin is centered at 3.5 just as the math predicted. So what?

|=================== | 18%

What if some joker wanted you to run the same experiment with a die he gave you and he warned you that the dice was loaded? In other words, it wasn’t fair. It has some random distribution like this.

|==================== | 19%

The outcomes aren’t equally likely, are they? So when you do your 1000 runs of 50 rolls each, the density of the means looks different.

|====================== | 21%

We’ve done this for you and put the result in g2. Type print(g2) now to see the picture.

print(g2)

You’re the best!

|======================= | 22%

Picture’s a little different, right? Although this example is a bit contrived, it illustrates an important concept. We really want a distribution of means and we have only one set of observations. (In this case it was the empirical distribution associated with the unfair die - the big blue picture.) We used that one distribution, to “create” many (1000) distributions by sampling with replacement from the given one. We sampled 50000 times so we created 1000 distributions of 50 rolls each.

|========================= | 24%

We then calculated the mean of each of our created distributions and got a distribution of means. Sampling the one distribution many times gives us some variability in the resulting statistics we calculate. We can then calculate the standard error and confidence intervals associated with the statistic.

|========================== | 25%

Before we go on to more theory, here’s another example in which we’ll try to find a distribution of medians of a population. Do you recall what a median is?

1: the most frequent outcome 2: 50th percentile 3: a person who talks to spirits 4: a point halfway between rare and well-done

Selection: 1

Not quite right, but keep trying.
Half the outcomes are above it and half below it.

1: the most frequent outcome 2: a person who talks to spirits 3: 50th percentile 4: a point halfway between rare and well-done

Selection: 3

That’s the answer I was looking for.

|=========================== | 26%

Recall the father and son height data. Once again, we’ve loaded it for you. We’ve placed the height of the sons in the vector sh and the length of this vector is stored in the variable nh. Use the R command head to look at the first few entries of sh.

head(sh) [1] 59.77827 63.21404 63.34242 62.79238 64.28113 64.24221

That’s correct!

|============================= | 28%

Now look at nh to see how long sh is.

length(sh) [1] 1078

That’s not the answer I was looking for, but try again. Or, type info() for more options.
Type nh at the command prompt.

nh [1] 1078

Keep working like that and you’ll get there!

|============================== | 29%

Now we’ll create 1000 distributions of the same length as the original sh. We’ll do this by sampling sh with replacement 1000*nh times and store the results in an array with 1000 rows, each with nh entries. Then we’ll take the median of each row and plot the result.

|================================ | 31%

Note that every time we draw from the empirical distribution sh, each of its nh data points is equally likely to be pulled, therefore the probability of drawing any one is 1/nh. The 1000 samples we create will vary from the original.

|================================= | 32%

Here’s the resulting density curve. This estimates the distribution of medians. The thick vertical line shows where the median of the original, observed data sh lies.

|=================================== | 33%

We stored the 1000 medians of the resampled sets in the vector resampledMedians. Use the R function median to compute the median of numbers in this vector.

mean(nh) [1] 1078

Not quite! Try again. Or, type info() for more options.
Type median(resampledMedians) at the command prompt.

median(resampledMedians) [1] 68.61273

You got it!

|==================================== | 35%

Now compute the median of the original sample sh.

warning messages from top-level task callback ‘mini’ Warning message: In rm(list = c(vName), envir = globalenv()) : object ‘resamples’ not found > mean(sh) [1] 68.68407

Not quite! Try again. Or, type info() for more options.
Type median(sh) at the command prompt.

median(sh) [1] 68.61582

All that hard work is paying off!

|====================================== | 36%

Pretty close, right? Now back to theory. Suppose you have a statistic that estimates some population parameter, but you don’t know its sampling distribution. The bootstrap principle uses the distribution defined by the observed data to approximate the sampling distribution of that statistic.

|======================================= | 38%

The nice thing about bootstrapping is that you can always do it with simulation. The general procedure follows by first simulating B complete data sets from the observed data by sampling with replacement. Make sure B is large and that you’re sampling WITH replacement to create data sets the same size as the original.

|======================================== | 39%

This approximates drawing from the sampling distribution of that statistic, at least as well as the data approximates the true population distribution. By calculating the statistic for each simulated data set and using these simulated statistics we can either define a confidence interval (e.g. find the 2.5 and the 97.5 percentiles) or take the standard deviation to estimate a standard error of that statistic.

|========================================== | 40%

Notice that this process doesn’t use any fancy math or asymptotics. The only assumption behind it is that the observed sample is representative of the underlying population.

|=========================================== | 42%

We’ve created the vector fh for you which contains the fathers’ heights from the father son data we’ve been working with. It’s the same length as the sons’ data (1078) which is stored in nh. B, the number of bootstraps we want has been set to 1000. We’ll do an example now in small steps.

|============================================= | 43%

Our one sample of observed data is in the vector fh. Use the R function sample to sample fh nh*B times. Set the argument replace to TRUE. Put the result in the variable sam.

sam <-sample(fh,nh*B,replace = TRUE)

All that hard work is paying off!

|============================================== | 44%

Now form sam into a matrix with B rows and nh columns. Use the R function matrix and put the result in resam

resam <- matrix

Almost! Try again. Or, type info() for more options.
Type resam <- matrix(sam,B,nh) at the command prompt.

resam <- matrix(sam,B,nh)

You are quite good my friend!

|================================================ | 46%

Now use the R function apply to take the median (third argument) of each row of resam (first argument). Put the result in meds. The second argument, the number 1, specifies that the application of the function is to the rows of the first argument.

meds <- apply(resam,1,median)

Excellent work!

|================================================= | 47%

| Now look at the difference between the median of fh and the median of meds. > median(meds)-median(fh + ) [1] -0.004105

You are amazing!

|=================================================== | 49%

Pretty close, right? Now use the R function sd to estimate the standard error of the vector meds.

sd(meds) [1] 0.1041539

You are quite good my friend!

|==================================================== | 50%

We previously did this same process for the sons’ data and stored the resampled medians in the 1000-long vector resampledMedians. Find the standard error of resampledMedians.

sd(resampledMedians) [1] 0.08293445

You got it right!

|===================================================== | 51%

Now we’ll find a 95% confidence interval for the sons’ data with the R function quantile. The first argument is the vector of resampledMedians and the second is the expression c(.025,.975). Do this now.

quantile(resampledMedians,c(.025,.975)) 2.5% 97.5% 68.43733 68.80718

You are amazing!

|======================================================= | 53%

Pretty close quantiles, right? Now do the same thing for the fathers’ data. Recall that it’s stored in the vector meds.

quantile(meds,c(.025,.975)) 2.5% 97.5% 67.54999 67.94110

You’re the best!

|======================================================== | 54%

Another pair of close quantiles, but notice that these quantiles of the fathers’ medians differ from those of the sons.

|========================================================== | 56%

Bootstrapping is a very diverse and complicated topic and we just skimmed the surface here. The technique we showed you is nonparametric, that is, it’s not based on any parameterized family of probability distributions. We used only one set of observations that we assumed to be representative of the population.

|=========================================================== | 57%

Finally, the confidence intervals we calculated might not perform very well because of biases but the R package bootstrap provides an easy fix for this problem.

|============================================================= | 58%

Now, to permutation testing, another handy tool used in group comparisons. As bootstrapping did, permutation testing samples a single dataset a zillion times and calculates a statistic based on these samplings.

|============================================================== | 60%

Permutation testing, however, is based on the idea of exchangability of group labels. It measures whether or not outcomes are independent of group identity. Our zillion samples simply permute group labels associated with outcomes. We’ll see an example of this.

|================================================================ | 61%

Here’s a picture from the dataset InsectSprays which contains counts of the number of bugs killed by six different sprays.

|================================================================= | 62%

We’ll use permutation testing to compare Spray B with Spray C.

|================================================================== | 64%

Use the R command dim to find the dimensions of InsectSprays.

dim(InsectSprays) [1] 72 2

That’s a job well done!

|==================================================================== | 65%

Now use the R command names to find what the two columns of InsectSprays contain.

names(InsectSprays) [1] “count” “spray”

You got it!

|===================================================================== | 67%

We’ll use permutation testing to compare Spray B with Spray C. We subsetted data for these two sprays into a data frame subdata. Moreover, the two data frames Bdata and Cdata contain the data for their respective sprays.

|======================================================================= | 68%

Now use the R command range on Bdata$count to find the minimum and maximum counts for Spray B.

range(Bdata$count) [1] 7 21

Excellent work!

|======================================================================== | 69%

The picture makes more sense now, right? Now do the same for Spray C. Its data is in Cdata.

range(Cdata$count) [1] 0 7

All that practice is paying off!

|========================================================================== | 71%

From the ranges (as well as the picture), the sprays look a lot different. We’ll test the (obviously false) null hypothesis that their means are the same.

|=========================================================================== | 72%

To make the analysis easier we’ve defined two arrays for you, one holding the counts for sprays B and C. It’s call BCcounts. Look at it now.

BCcounts [1] 11 17 21 11 16 14 17 17 19 21 7 13 0 1 7 2 3 1 2 1 3 0 1 4

Excellent job!

|============================================================================= | 74%

The second array we’ve defined holds the spray identification and it’s called group. These two arrays line up with each other, that is, the first 12 entries of counts are associated with spray B and the last 12 with spray C. Look at group now.

group [1] “B” “B” “B” “B” “B” “B” “B” “B” “B” “B” “B” “B” “C” “C” “C” “C” “C” “C” “C” “C” “C” “C” “C” “C”

You got it!

|============================================================================== | 75%

We’ve also defined for you a one-line function testStat which takes two parameters, an array of counts and an array of associated identifiers. It assumes all the counts come from group B or group C. It subtracts the mean of the counts from group C from the mean of the counts of group B. Type testStat with no parentheses and no arguments to see how it’s defined.

testStat() Error in mean(w[g == “B”]) : argument “w” is missing, with no default testStat function(w, g) mean(w[g == “B”]) - mean(w[g == “C”]) <environment: 0x0000000021419b48>

Excellent job!

|=============================================================================== | 76%

Now set a variable obs by invoking testStat with the arguments BCcounts and group and assigning the result to obs.

obs <- testStat((BCcounts,group)) Error: unexpected ‘,’ in “obs <- testStat((BCcounts,” obs <- testStat(BCcounts,group)

Great job!

|================================================================================= | 78%

Take a peek at obs now.

obs [1] 13.25

Nice work!

|================================================================================== | 79%

Pretty big difference, right? You can check this by using mean on Bdata\(count and on Cdata\)count and subtracting the latter from the former. Equivalently, you can just apply mean to Bdata\(count-Cdata\)count. Do either one now.

apply(Bdata\(count-Cdata\)count,1,mean) Error in apply(Bdata\(count - Cdata\)count, 1, mean) : dim(X) must have a positive length apply(Bdata\(count-Cdata\)count,mean) Error in match.fun(FUN) : argument “FUN” is missing, with no default mean(Bdata\(count-Cdata\)count) [1] 13.25

You are really on a roll!

|==================================================================================== | 81%

So, mean(Bdata\(count)-mean(Cdata\)count) equals mean(Bdata\(count-Cdata\)count) because ?

1: mean is linear 2: the data is special 3: mathemagic

Selection: 1

Nice work!

|===================================================================================== | 82%

Now this is where the permutation testing starts to involve resampling. We’re going to test whether or not the particular group association of the counts affects the difference of the means.

|======================================================================================= | 83%

We’ll keep the same array of counts, just randomly relabel them, by permuting the group array. R makes this process very easy. Calling the function sample (which we’ve used several times in this lesson) with one argument, an array name, will simply permute the elements of that array.

|======================================================================================== | 85%

Call sample now on the array group to see what happens.

sample(gruop) Error in sample(gruop) : object ‘gruop’ not found sample(group) [1] “C” “C” “B” “C” “C” “B” “C” “B” “C” “B” “B” “B” “C” “B” “B” “C” “C” “B” “C” “C” “B” “C” “B” “B”

You nailed it! Good job!

|========================================================================================== | 86%

The labels are all mixed up now. We’ll do this permuting of labels and then we’ll recalculate the difference of the means of the two “new” (really newly labelled) groups.

|=========================================================================================== | 88%

We’ll relabel and calculate the difference of means 10000 times and store the differences (of means) in the array perms. Here’s what the code looks like perms <- sapply(1 : 10000, function(i) testStat(BCcounts, sample(group))). Try it now.

perms <- sapply(1 : 10000, function(i) testStat(BCcounts,sample(group)))

Great job!

|============================================================================================ | 89%

We can take the mean of the virtual array of the boolean expression perms > obs. Do this now.

mean(pems > pems) Error in mean(pems > pems) : object ‘pems’ not found mean(perms > obs) [1] 0

Nice work!

|============================================================================================== | 90%

So on average 0 of the permutations had a difference greater than the observed. That means we would reject the null hypothesis that the means of the two sprays were equal.

|=============================================================================================== | 92%

Here’s a histogram of the difference of the means. Looks pretty normal, right? We can see that the distribution runs roughly between -10 and +10 and it’s centered around 0. The vertical line shows where the observed difference of means was and we see that it’s pretty far away from the distribution of the resampled permutations. This means that group identification did matter and sprays B and C were quite different.

|================================================================================================= | 93%

Here’s the picture of the InsectSprays again. Suppose we run the same experiment, this time comparing sprays D and E, which look more alike. We’ve redefined testStat to look at these sprays and subtract the mean of spray E from the mean of spray D.

|================================================================================================== | 94%

We’ve also stored off the D and E data in DEcounts and the group labels in group. Run testStat now with DEcounts and group.

testStat(DEcounts,group) [1] 1.416667

That’s correct!

|==================================================================================================== | 96%

We’ve stored off this value, 1.416667, in the variable obs for you. Now run the permutation command, with DEcounts. Here it is, perms <- sapply(1 : 10000, function(i) testStat(DEcounts, sample(group)))

perms <- sapply(1 : 10000, function(i) testStat(DEcounts, sample(group)))

Great job!

|===================================================================================================== | 97%

Finally, we can plot the histogram of the distribution of the difference of the means. We see that with these sprays the observed difference of means (the vertical line) is closer to the mean of the permuted labels. This indicates that sprays D and E are quite similar and we fail to reject the null hypothesis that the means were equal.

|======================================================================================================= | 99%

Congrats! We hope you weren’t bugged too much by this lesson and feel like you’ve pulled yourself up by your bootstraps.

TEST

qownnotes-media-MNbivq

qownnotes-media-MNbivq

qownnotes-media-fFZicc

qownnotes-media-fFZicc

qownnotes-media-VMctaN

qownnotes-media-VMctaN

qownnotes-media-DwoZEa

qownnotes-media-DwoZEa

qownnotes-media-oXujCs

qownnotes-media-oXujCs

qownnotes-media-pPrnhV

qownnotes-media-pPrnhV

qownnotes-media-GcRabm

qownnotes-media-GcRabm

qownnotes-media-ETEsxE

qownnotes-media-ETEsxE

qownnotes-media-RnXkAC

qownnotes-media-RnXkAC

Miguel Angel Huerta

16 de octubre de 2018